Theory Digest — Apr 21, 2026

1. Motivation (2–3 sentences): Understanding generalization in neural networks trained at large learning rates (“edge of stability”) where optimization exhibits chaotic behavior is a central open problem in machine learning. Classical convex optimization theory breaks down when the largest Hessian eigenvalue exceeds the stability threshold 2/η, yet empirically this regime often yields improved generalization.

2. Mathematical setup: Consider supervised learning with parameter vector $w \in \mathbb{R}^d$ and population risk:

\[R(w) := \mathbb{E}_{Z\sim\mu}[\ell(w, Z)]\]

where $\mu$ is the data distribution and $\ell: \mathbb{R}^d \times \mathcal{Z} \rightarrow \mathbb{R}$ is the loss function. Training minimizes empirical risk:

\[\hat{R}_S(w) := \frac{1}{n}\sum_{i=1}^n \ell(w, Z_i)\]

over dataset $S = {Z_i}_{i=1}^n \sim \mu^{\otimes n}$ using stochastic optimization. The optimization dynamics are modeled as a random dynamical system (RDS) with cocycle $\phi: \mathbb{N}_0 \times \Omega \times \mathbb{R}^d \rightarrow \mathbb{R}^d$ over metric dynamical system $(\Omega, \mathcal{F}, P, \theta)$.

Assumptions:

1. Loss function is bounded: $0 \leq \ell(w,z) \leq B$
2. Loss is L-Lipschitz continuous in $w$
3. RDS admits unique compact random pullback attractor $A_S(\omega)$

4. Non-singularity and integrability conditions on the Jacobian $D\phi(1,\omega,x)$

   3. Toy example: Consider gradient descent with fixed learning rate $\eta$ on a quadratic loss where the Hessian has eigenvalues $\alpha_1 > 2/\eta$ and $\alpha_2 < 2/\eta$. The update Jacobian $I - \eta \nabla^2 \hat{R}_S(w)$ has singular values $

   1-\eta\alpha_i

   $, giving expansion along the first eigendirection when $\alpha_1 > 2/\eta$.

   4. Formal objective: The goal is to bound the worst-case generalization gap over the random attractor:

   \[G(A_S(\omega)) := \sup_{w \in A_S(\omega)} R(w) - \hat{R}_S(w)\]

The key innovation is the Sharpness Dimension (SD), defined via RDS sharpness measures.

RDS Sharpness of Order k: For RDS Jacobian $D\phi(1,\omega,w)$ with singular values $\sigma_1(\omega,w) \geq \cdots \geq \sigma_d(\omega,w)$:

\[\lambda_k := \mathbb{E}\left[\sup_{w \in A(\omega)} \ln \sigma_k(\omega,w)\right]\]

Sharpness Dimension: Given ordered sharpness indices $\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_d$, define:

\[j^* := \max\left\{i \in \{1,\ldots,d\} : \sum_{k=1}^i \lambda_k \geq 0\right\}\]

Then:

\[\dim_S A := \begin{cases}\] \[j^* + \frac{\sum_{i=1}^{j^*} \lambda_i}{|\lambda_{j^*+1}|} & \text{if } 1 \leq j^* < d \\\] \[d & \text{if } j^* = d \\\] \[0 & \text{if } \lambda_1 < 0\] \[\end{cases}\]

Algorithm steps:

1. Train network using stochastic optimizer (e.g. SGD)
2. After convergence, estimate Hessian spectrum via stochastic Lanczos quadrature (SLQ)
3. Compute RDS sharpness $\lambda_k$ from Jacobian $I - \eta \nabla^2 \hat{R}_S(w)$ singular values

4. Calculate SD using the formula above

   Applied to toy example: For the quadratic case with eigenvalues $\alpha_1 > 2/\eta > \alpha_2$, we get $\lambda_1 = \ln

   1-\eta\alpha_1

   > 0$ and $\lambda_2 = \ln

   1-\eta\alpha_2

   < 0$. Since $\lambda_1 > 0 > \lambda_1 + \lambda_2$, we have $j^* = 1$ and:

   \[\dim_S A = 1 + \frac{\lambda_1}{|\lambda_2|} < d = 2\]

Step 1 — Prior work:

• Cohen et al. (2021) “Gradient descent on neural networks typically occurs at the edge of stability” - introduced the edge of stability phenomenon where largest Hessian eigenvalue exceeds 2/η
• Camuto et al. (2021) “Fractal structure and generalization properties of stochastic optimization algorithms” - linked fractal geometry to generalization but only under contraction, not applicable at EoS
• Simsekli et al. (2020) “Hausdorff dimension, heavy tails, and generalization in neural networks” - established topology-based generalization bounds via trajectory analysis

Step 2 — Delta: This paper extends generalization theory to the unstable/chaotic regime by:

1. modeling optimization as random dynamical systems with pullback attractors instead of point convergence
2. introducing Sharpness Dimension as a Lyapunov-inspired complexity measure capturing effective dimensionality of expanding directions

3. proving generalization bounds in terms of attractor geometry rather than individual solutions.

   Step 3 — Theory-specific assessment:

   • The main theorem connecting SD to generalization gap is technically sound but follows expected lines - it reduces to covering number arguments via Minkowski dimension bounds
   • The proof technique of using ellipsoidal coverings determined by the Hessian spectrum is relatively standard in fractal geometry
   • No lower bounds are established for the Sharpness Dimension bound, so tightness is unknown
   • The RDS formulation, while mathematically elegant, doesn’t fundamentally change the analysis compared to existing trajectory-based approaches

   Verdict: INCREMENTAL — Solid theoretical extension of existing fractal dimension generalization bounds to the edge of stability regime, but the core techniques and results are predictable refinements of prior work.

The main proof strategy consists of three stages:

1. Dimension comparison: Show that Sharpness Dimension upper bounds Minkowski dimension of the attractor:

   \[\dim_M A_S(\omega) \leq \dim_S A_S\]

   This is achieved by constructing ellipsoidal coverings of $A_S(\omega)$ where the ellipsoid axes are determined by the RDS sharpness values $\lambda_1, \ldots, \lambda_d$.

2. Volume growth control: The key inequality leverages the bounded distortion assumption (Assumption 4.4) to control how $j$-dimensional volumes evolve under the dynamics:

   \[\mathbb{E}\left[\sup_{x \in A(\omega)} \ln \|D\phi(m,\omega,x)\|_j - \inf_{x \in A(\omega)} \ln \|D\phi(m,\omega,x)\|_j\right] = o(m)\]

   where $|A|_j := \sigma_1(A) \cdots \sigma_j(A)$ is the maximal $j$-volume expansion factor.

3. Covering number bound: Using the ellipsoidal covering construction and multiplicative ergodic theorem, the authors bound the $\epsilon$-covering number:

   \[\mathcal{N}(A_S(\omega), \epsilon) \leq C \epsilon^{-\dim_S A_S}\]

   for some constant $C > 0$.

4. Generalization bound: The final step applies existing PAC-Bayesian results (Dupuis et al. 2024) connecting covering numbers to generalization via:

   \[G(A_S(\omega)) \leq 2L\delta + 2B\sqrt{\frac{4 \dim_S A_S \log(1/\delta)}{n}} + \text{lower order terms}\]

   The proof relies heavily on the technical machinery of random dynamical systems theory, particularly the multiplicative ergodic theorem and properties of random pullback attractors.

Datasets: MNIST (3-layer MLP with 12,960 parameters, 5-layer MLP with 278,800 parameters), WikiText-2 (GPT-2 with 124M parameters), arithmetic modulo 97 (2-layer MLP, grokking experiments).

Baselines: Persistent homology dimensions (PH-Dim), α-weighted lifetime sums (E_α), Hessian trace, classical sharpness (largest Hessian eigenvalue), trajectory-based topological measures.

Key results:

• Small MLPs: SD achieves correlation 0.83 with generalization gap vs. 0.66 for best baseline
• Larger networks: SD-SLQ maintains strong correlations (τ = 0.4) while exact methods become intractable
• GPT-2: SD variants show most consistent positive correlations across SGD, SGD+momentum, AdamW
• Grokking: SD captures phase transitions during sudden generalization with correlations 0.69-0.88

Implementation: Uses stochastic Lanczos quadrature (SLQ) with 500 runs of 100 steps each for large-scale Hessian spectrum estimation. Computational complexity scales as O(md log d) per iteration where m is number of Lanczos steps.

Limitations:

1. TECHNICAL: Requires estimation of full Hessian spectrum, computationally demanding for very large models despite SLQ approximations
2. TECHNICAL: Bounded distortion assumption (Assumption 4.4) may be difficult to verify in practice and could be restrictive for some architectures
3. NATURAL: Theory currently limited to discrete-time dynamics; extension to continuous-time stochastic differential equations would be valuable
4. RESTRICTIVE: Non-singularity condition requiring inf σ_d(Dφ) > 0 excludes degenerate cases and may not hold for all network architectures

5. TECHNICAL: Integrability assumptions on second derivatives needed for C² regularity may be violated in practice

   Open problems:

6. Extend theory to adaptive optimizers like Adam/AdamW which exhibit different spectral properties than SGD
7. Establish lower bounds for the Sharpness Dimension generalization bound to determine tightness

1. Motivation (2–3 sentences): In Bayesian optimization with Gaussian process surrogates, Thompson sampling must evaluate function samples at discrete candidate points due to computational intractability. High-dimensional problems suffer from the curse of dimensionality where candidate points become exponentially sparse, severely degrading the quality of Thompson sampling approximations.

2. Mathematical setup: Consider optimizing $F(x): \mathbb{R}^d \to \mathbb{R}$ over compact domain $X \subset \mathbb{R}^d$ with noisy observations $y = F(x) + \sigma_n \epsilon$ where $\epsilon \sim \mathcal{N}(0,1)$. The GP posterior after $n_t$ observations $D_t = {(x_i, y_i)}_{i=1}^{n_t}$ is:

   \[p(f | D_t) = \text{GP}(\mu_t(x), k_t(x,x'))\]

   where:

   \[\mu_t(x) = k(x, X_t)^T K^{-1} y\] \[k_t(x,x') = k(x,x') - k(x,X_t)^T K^{-1} k(X_t, x')\]

   Assumptions:

   1. GP prior $p(f) = \text{GP}(0, k(x,x’))$ with once-differentiable kernel $k$
   2. Kernel hyperparameters estimated via type-II MLE
   3. Compact search domain $X \subset \mathbb{R}^d$
   4. Bounded noise variance $\sigma_n^2$

   Thompson sampling draws $f \sim p(f

   D_t)$ and selects $x_{t+1} = \arg\max_{x \in X} f(x)$, but this requires discretization over candidate points $\tilde{X} = {x_1^*, \ldots, x_M^*}$.

3. Toy example: When $d=2$, $X=[0,1]^2$, and using $M=100$ Sobol candidate points, the candidate density is only $100/1 = 100$ points per unit area. For $d=10$ with same budget, density drops to $100/1 = 100$ points per unit volume, but the volume is still 1, so effective density becomes critically sparse for finding function maxima.

4. Formal objective: Find the next evaluation point by solving:

   \[x_{t+1} = \arg\max_{x^* \in \tilde{X}} f(x^*)\]

   where $f \sim p(f

   D_t)$ and $\tilde{X}$ is adaptively constructed to maximize the quality of this discrete optimization.

ACTS generates candidate points within a gradient-aligned subspace to increase density in promising regions. The algorithm:

1. Sample incumbent: $x_0 = \arg\max_{x \in D_t} \mathbb{E}[f(x)

   D_t]$

2. Sample gradient: $\nabla f(x_0) \sim p(\nabla f(x_0)

   D_t)$ using joint GP model:

   \[\begin{pmatrix} y_t \\ \nabla f(x_0) \end{pmatrix} \sim \mathcal{N}\left(0, \begin{pmatrix} K + \sigma_n^2 I & k(X_t, x_0)_\nabla^T \\ \nabla k(x_0, X_t) & \nabla k(x_0, x_0)_\nabla^T \end{pmatrix}\right)\]

3. Construct adaptive search space as axis-aligned cone:

   \[T_{\nabla f(x_0)} = \{x_0 + v \odot \nabla f(x_0) \mid 0 \preceq v \in \mathbb{R}^d\} \cap X\]

4. Generate candidates: $\tilde{X}_t \sim \pi(\tilde{X}

   \nabla f(x_0))$ by applying base policy $\pi_0$ to $T_{\nabla f(x_0)}$

5. Sample function values: $f_{\tilde{X}_t} \sim p(f_{\tilde{X}_t}

   \nabla f(x_0), D_t)$

6. Select maximizer: $x_{t+1} = x^*_{i_{\max}}$ where $i_{\max} = \arg\max_i [f_{\tilde{X}_t}]_i$

   For the RAASP base policy, ACTS modifies the dimension selection probabilities:

   \[[b_i]_j \sim \text{Bernoulli}\left(\min\left\{\frac{20[\nabla f(x_0)]_j^2}{||\nabla f(x_0)||_2^2}, 1\right\}\right)\]

   Applied to toy example: With $d=2$, $x_0 = (0.5, 0.5)$, and sampled gradient $\nabla f(x_0) = (0.3, -0.1)$, the search space becomes $T = {(0.5, 0.5) + (v_1 \cdot 0.3, v_2 \cdot (-0.1)) \mid v_1, v_2 \geq 0} \cap [0,1]^2$, which has volume approximately $0.5 \times 0.5 \times 0.3 \times 0.1 = 0.0075$ compared to the full domain volume of 1, yielding $100/0.0075 \approx 13,333$ effective candidate density.

Step 1 — Prior work:

• “Scalable global optimization via local Bayesian optimization (TuRBO)” (Eriksson et al., 2019): Uses trust regions and RAASP candidate policy for high-dimensional BO, achieving strong performance by restricting search to axis-aligned subspaces.
• “Cylindrical Thompson sampling for high-dimensional Bayesian optimization” (Rashidi et al., 2024): Generates candidates in spherical regions around incumbents to improve local exploration in Thompson sampling.
• “Efficiently sampling functions from Gaussian process posteriors” (Wilson et al., 2020): Pathwise sampling approach avoiding candidate points by approximating GP with finite basis functions.

Step 2 — Delta: This paper introduces adaptive candidate point construction that couples the candidate set selection with the posterior sample itself. Key innovations:

1. Joint sampling of function values and gradients to inform candidate placement
2. Gradient-aligned cone construction that dramatically reduces search volume while maintaining global consistency

3. Integration with existing candidate policies like RAASP through modified dimension selection.

   Step 3 — Theory-specific assessment:

   • The main result (global consistency, Theorem 1) is not surprising given that the sampled gradients retain randomness and don’t fully collapse posterior uncertainty. The proof follows standard concentration arguments.
   • The proof technique is routine, applying standard GP concentration bounds and compactness arguments. No genuinely novel technical insights.
   • No tightness analysis provided. The volume reduction (factor of $2^{100} \approx 10^{30}$ in 100D) is impressive but the connection to optimization performance bounds is not rigorously established.

   The adaptive candidate construction is conceptually neat but represents an incremental improvement over existing sparse/local strategies rather than a fundamental breakthrough.

   Verdict: INCREMENTAL — solid engineering contribution with reasonable theoretical backing, but the core insight (align candidates with gradients) is fairly natural and the theoretical analysis is standard.

The main theoretical result (Theorem 1) proves global consistency using contradiction and concentration arguments:

1. Proof by contradiction: Assume there exists $\epsilon > 0$ such that ACTS never queries within $\epsilon$ of the global optimum $x^*$.

2. Key concentration bound: For any unobserved point $x_u$, the posterior variance is lower bounded:

   \[\sigma_t^2(x_u) \geq \delta > 0\]

   for some $\delta$ depending on kernel properties and the separation from observed points.

3. Gradient non-degeneracy: Since posterior covariance never fully collapses, sampled gradients satisfy:

   \[P(||\nabla f(x_0)||_2 > 0) = 1\]

   almost surely, ensuring the search space $T_{\nabla f(x_0)}$ remains well-defined.

4. Coverage argument: Under assumption A6, the candidate distribution has full support over $T_{\nabla f(x_0)} \cap X_\epsilon$, where $X_\epsilon$ is an $\epsilon$-covering of $X$.

5. Infinite occurrence: The event that the adaptive search space contains points near $x^*$ occurs infinitely often with probability 1, leading to eventual exploration.

   The proof is technically sound but follows standard templates for GP-based optimization convergence. The main technical lemma bounds posterior variance decay, which is routine given bounded kernels and compact domains. No novel concentration inequalities or sophisticated proof techniques are introduced.

Datasets: Medium-dimensional robotics tasks (Lunar Lander 12D, Robot Pushing 14D, Swimmer 16D, Hopper 32D), high-dimensional benchmarks (Rover 60D, MOPTA08 124D, SVM hyperparameters 388D, LassoBench 180-1000D, GuacaMol molecular design 256D).

Baselines: Thompson sampling variants (RAASP, Cylindrical TS, Pathwise TS), state-of-the-art methods (SAASBO, BAxUS, LogEI), with/without TuRBO trust regions.

Key numbers: ACTS achieves top performance (within 2 standard errors) on most benchmarks. On MOPTA08, ACTS reaches -215.81 vs RAASP’s -225.12. Search space volume reduction by factors up to $10^{30}$ in high dimensions. Computational overhead is minimal with $O(M^3)$ complexity maintained.

Batch optimization with $q=10,50,100$ shows consistent performance. Wallclock timing shows ACTS adds only 20% overhead over RAASP (0.53s vs 0.44s per iteration on Rover 60D).

Limitations:

1. Requires once-differentiable kernel - NATURAL (standard requirement, satisfied by RBF, Matérn kernels)
2. Assumes gradient-aligned regions contain local maxima - TECHNICAL (reasonable for smooth functions but may fail with very small lengthscales)
3. Axis-aligned cone geometry is somewhat ad-hoc - TECHNICAL (authors test alternatives in appendix, current choice works well empirically)
4. Single-step gradient adaptation may be myopic - TECHNICAL (could extend to multi-step but adds computational cost)

5. Limited to continuous domains - NATURAL (standard BO assumption)

   Open problems:

6. Theoretical characterization of when gradient-aligned search spaces contain near-optimal regions, potentially connecting to local smoothness assumptions and GP lengthscale parameters.
7. Extension to multi-step gradient following (mimicking gradient ascent on GP samples) while maintaining computational efficiency and theoretical guarantees.

1. Motivation: This work extends Nesterov’s accelerated gradient descent beyond smooth strongly convex optimization to the broader class of strongly monotone, Lipschitz operators. This matters for minimax optimization problems and monotone inclusion problems that arise in game theory, robust optimization, and variational inequalities.

2. Mathematical setup: Consider the monotone inclusion problem

   \[\text{find } z \in \mathbb{R}^d \text{ such that } \mathbb{T}(z) = 0,\]

   where $\mathbb{T}: \mathbb{R}^d \to \mathbb{R}^d$ is maximal, strongly monotone, and Lipschitz. Assumptions:

   1. $\mathbb{T}$ is $\mu$-strongly monotone: $\langle w - w’, z - z’ \rangle \geq \mu |z - z’|^2$ for all $z, z’ \in \mathbb{R}^d$, $w \in \mathbb{T}(z)$, $w’ \in \mathbb{T}(z’)$
   2. $\mathbb{T}$ is $L$-Lipschitz: $|w - w’| \leq L |z - z’|$
   3. $\mathbb{T}$ admits Asplund decomposition $\mathbb{T} = \partial\phi + \mathbb{S}$ where $\phi$ is convex and $\mathbb{S}$ is acyclic monotone

3. Toy example: When $d = 2$ and $L(x,y) = x^2 - y^2 + \sin x \sin y$, the saddle operator is

   \[\mathbb{T}(x,y) = \begin{pmatrix} 2x + \cos x \sin y \\ 2y - \sin x \cos y \end{pmatrix}.\]

   This decomposes as $\mathbb{T} = \nabla\phi + \mathbb{S}$ with non-bilinear coupling, illustrating the core challenge of handling general monotone structure beyond gradients.

4. Formal objective: Find $z^*$ such that $\mathbb{T}(z^*) = 0$ with iteration complexity

   \[O\left(\sqrt{\frac{L_\phi}{\mu} + \frac{L_\mathbb{S}^2}{\mu^2}} \log \frac{1}{\epsilon}\right).\]

The Nesterov acceleration with Operator Decomposition (NOD) algorithm operates on the decomposition $\mathbb{T} = \nabla\phi + \mathbb{S}$:

1. Gradient step: $z_{k+1} = \tilde{z}_k - \eta(\nabla\phi(\tilde{z}_k) + \mathbb{S}(\hat{z}_k))$
2. Standard momentum: $\tilde{z}_{k+1} = z_{k+1} + \tau(z_{k+1} - z_k)$ where $\tau = \frac{1-\sqrt{\eta\mu}}{1+\sqrt{\eta\mu}}$

3. Extended momentum: $\hat{z}_{k+1} = \tilde{z}_{k+1} + \theta(\tilde{z}_{k+1} - \tilde{z}_k)$ where $\theta = \frac{1}{\sqrt{\eta\mu}} - 1$

   The key insight is using different extrapolation points: $\tilde{z}_k$ for the cyclically monotone component $\nabla\phi$ and $\hat{z}_k$ (with larger extrapolation) for the acyclic monotone component $\mathbb{S}$.

   Applied to toy example: For $L(x,y) = x^2 - y^2 + \sin x \sin y$, NOD evaluates $\nabla\phi$ at the standard Nesterov extrapolation point $\tilde{z}_k$, but evaluates the coupling term $\mathbb{S}$ at the more aggressively extrapolated point $\hat{z}_k$.

Step 1 — Prior work:

• “A unified analysis of extra-gradient and optimistic gradient methods for saddle point problems” (Mertikopoulos et al., 2019): achieved $O(\sqrt{L/\mu} \log(1/\epsilon))$ rates for bilinear coupling only
• “Accelerated gradient methods for nonconvex-concave and convex-nonconcave minimax problems” (Lin et al., 2020): similar rates but restricted to specific coupling structures
• “On the convergence rate of training recurrent neural networks” (Chen et al., 2018): developed continuous-time analysis for accelerated methods

Step 2 — Delta: This paper extends acceleration to general strongly monotone operators via Asplund decomposition, moving beyond bilinear coupling to arbitrary monotone structure. The algorithm uses decomposition-aware extrapolation with different momentum parameters for cyclically monotone vs acyclic components.

Step 3 — Theory-specific assessment:

• Main theorem is somewhat predictable: rates match existing bilinear results when specialized
• Proof technique is routine: assembles known concentration and smoothness arguments with careful bookkeeping for dual extrapolation
• Bounds appear to match structure of the problem but no lower bounds are established to assess tightness

The Asplund decomposition characterization for $L(x,y) = x^2 - y^2 + \sin x \sin y$ provides the first non-bilinear example, but this is more of a technical computation than conceptual breakthrough.

Verdict: INCREMENTAL — solid extension of known acceleration techniques to broader operator class using standard proof methods.

The proof uses a discrete Lyapunov function approach with careful decomposition-aware analysis:

1. Define modified Lyapunov function:

   \[\Psi_k = \mu\left\|\tilde{z}_k + \frac{1}{\sqrt{\eta\mu}}(\tilde{z}_k - z_k) - z^*\right\|^2 + 2\hat{\phi}(\tilde{z}_{k-1}) - \eta\|\nabla\phi(\tilde{z}_{k-1}) + \mathbb{S}(\hat{z}_{k-1})\|^2 + \text{correction terms}\]

2. Key contraction inequality through strong convexity:

   \[-2\langle\tilde{z}_k - z^*, \nabla\hat{\phi}(\tilde{z}_k)\rangle \leq -2\hat{\phi}(\tilde{z}_k) - \mu\|\tilde{z}_k - z^*\|^2\]

3. Monotonicity exploitation for acyclic component:

   \[-2\sqrt{\eta\mu}\langle\mathbb{S}(\hat{z}_k), \hat{z}_k - z^*\rangle \leq 0\]

4. Technical bound on extrapolation difference using Lipschitz property:

   \[\|\mathbb{S}(\hat{z}_{k-1}) - \mathbb{S}(\hat{z}_k)\| \leq L_\mathbb{S}\|\hat{z}_{k-1} - \hat{z}_k\|\]

5. Careful completion of squares and step size restriction:

   \[\eta \leq C \cdot \min\left\{\frac{1}{L_\phi}, \frac{\mu}{L_\mathbb{S}^2}\right\}\]

   where $C \approx 0.026118$ is a technical constant.

   The proof assembles standard convex analysis tools but requires intricate tracking of multiple momentum sequences and their interactions.

Purely theoretical. Empirical validation would require:

1. implementing Asplund decomposition computation for specific operator classes
2. comparing NOD against standard methods on strongly-convex-strongly-concave problems with non-bilinear coupling
3. demonstrating the $O(\sqrt{L_\phi/\mu + L_\mathbb{S}^2/\mu^2} \log(1/\epsilon))$ rate empirically.

Limitations:

1. Asplund decomposition is non-constructive via Zorn’s lemma — RESTRICTIVE (severely limits practical applicability)
2. Step size requires knowledge of both $L_\phi$ and $L_\mathbb{S}$ — TECHNICAL (standard in optimization literature)
3. Analysis requires maximal monotonicity — NATURAL (standard assumption for monotone inclusion problems)

4. Constants $B$ and $C$ appear suboptimal — TECHNICAL (likely artifact of proof technique)

   Open problems:

5. Develop constructive algorithms for computing Asplund decomposition beyond bilinear cases
6. Establish matching lower bounds for strongly monotone inclusion problems to assess optimality of the $O(\sqrt{L_\phi/\mu + L_\mathbb{S}^2/\mu^2})$ dependence