Adam's epsilon Is a Convergence-Critical Hyperparameter: Provable Divergence at Default Settings and the eps/Scale Phase Boundary

1. Introduction

Adam (Kingma and Ba, 2015) is the default optimizer of modern machine learning, yet its behaviour under the hyperparameters people actually ship is still not fully characterised. Reddi, Kale and Kumar (2018) showed that Adam can fail to converge on convex online problems, and a substantial literature has since mapped when this happens in the $({\beta }_1,{\beta }_2)$ plane (Zhang et al., 2022; 2026), produced positive rates under bounded gradients (Défossez et al., 2022), modified the update to remove the failure (Reddi et al., 2018; Taniguchi et al., 2024), or identified the biased limit points Adam does reach (Dereich and Jentzen, 2024).

One axis of the hyperparameter space has largely escaped this program: the stability constant $\epsilon$ that every implementation adds to the denominator, with default $10^{-8}$. Divergence constructions take $\epsilon =0$ or "sufficiently small" and treat this as harmless; positive results either require $\epsilon$ large (Zaheer et al., 2018) or degrade as $\epsilon \to 0$, e.g. by admitting only step sizes proportional to $\epsilon$. Practitioners meanwhile treat $\epsilon$ as a numerical fudge factor — while quietly tuning it by orders of magnitude in reinforcement learning and large-model training, and occasionally discovering that it changes outcomes qualitatively (Everett et al., 2024).

This paper closes that axis on the canonical counterexample family. We pose a question we have not found resolved anywhere:

Question. Does Adam with its exact shipped defaults $({\beta }_1,{\beta }_2,\epsilon )=(0.9,0.999,10^{-8})$, the positive $\epsilon$ included, converge on every convex stochastic problem with uniformly bounded stochastic gradients — under constant or $1/\sqrt{t}$ step sizes? If not, can the failure be repaired by tuning $\epsilon$ alone, and at what threshold?

Answer. No — and the threshold is a scale-coupled phase boundary. There is an explicit one-dimensional convex problem with bounded i.i.d. gradients on which default Adam ascends the objective almost surely at a linear rate, for every constant step size, and at a $\sqrt{T}$ rate under ${\alpha }_t=\alpha /\sqrt{t}$. The failure is governed by a single scalar drift $D({\beta }_1,{\beta }_2,\epsilon /\lambda )$, where $\lambda$ is the gradient scale: below a problem-dependent critical ratio $\epsilon /\lambda$ the drift is negative and Adam provably ascends; above it the drift is positive and Adam descends. Tuning $\epsilon$ alone rescues the entire family, uniformly over $({\beta }_1,{\beta }_2)$ — but no fixed $\epsilon$ is safe across gradient scales, so rescaling losses or rewards silently moves Adam across the boundary.

Contributions

C1 — exact reduction. On the canonical family, Adam's stationary moment estimates have the closed form $m=(1-{\beta }_1)(C+1)W_1-\lambda$ and $v=(1-{\beta }_2)(C^2-1)W_2+{\lambda }^2$, where $W_i={\sum }_{k\in S}{\beta }_i^k$ are geometric sums over the random lags $S$ of rare large gradients. All asymptotics reduce to the sign of $D=\mathbb{E}[m/(\sqrt{v}+\epsilon )]$ (Lemma 1).

C2 — almost-sure dichotomy. For every constant $\alpha >0$: $x_T/T\to -\alpha D$ almost surely, so $D<0$ means Adam ascends a convex objective at a linear rate — a qualitatively stronger failure than a regret lower bound or a noise floor: the iterate moves in the wrong direction forever, and $F(x_T)\to +\infty$. Under ${\alpha }_t=\alpha /\sqrt{t}$ the ascent persists at rate $\sqrt{T}$ (Theorem 2). A projected variant on $[-1,1]$ parks at the anti-minimizer up to $O(\alpha )$ (Proposition 3).

C3 — two-sided phase bounds; defaults provably diverge. Hand-checkable bounds: at the exact defaults, $D<0$ for an explicit instance ($C=10^4$, $p\approx 1.1\times 10^{-4}$) for every $\epsilon \le 17.4\lambda$ — covering $\epsilon =10^{-8}$ by nine orders of magnitude (Theorem 4, Corollary 5). Conversely $D>0$ whenever $\epsilon \ge \lambda [(C-1)/\mu -1]$, for all ${\beta }_1,{\beta }_2$ simultaneously (Theorem 4(ii)) — an $\epsilon$-analogue of the ${\beta }_2$-rescue of Zhang et al. (2022): on this family Adam is repaired by tuning $\epsilon$ alone, at fixed betas.

C4 — scale law. Adam run on the problem scaled by $\lambda$ is step-for-step identical to Adam run on the unit-scale problem with $\epsilon$ replaced by $\epsilon /\lambda$ (Lemma 1(ii)). Multiplying a loss by a constant — naive mixed-precision loss scaling, reward scaling, objective re-normalisation — can therefore flip Adam between convergence and divergence with every optimizer hyperparameter held fixed (Corollary 6). Adam-with-$\epsilon$ is not scale invariant, and the violation changes the convergence class, not just the constants.

C5 — certified, reproducible experiments. A dependency-free numpy implementation (Appendix B; all seeds fixed) that (a) evaluates $D$ by exact sampling of $(W_1,W_2)$ with rigorous truncation brackets plus Hoeffding confidence $1-10^{-6}$, certifying the sign at the defaults; (b) traces the empirical boundary ${\epsilon }^{\ast }(C)$ for three $({\beta }_1,{\beta }_2)$ regimes including the LLM-standard ${\beta }_2=0.95$; (c) cross-validates the theory by matching measured slopes of an exact bias-corrected Adam to the predicted $-\alpha D$; and (d) demonstrates the loss-scale flip with identical optimizer settings.

2. Related work and precise positioning

Divergence counterexamples. Reddi et al. (2018) construct online convex problems where Adam has non-vanishing regret; their argument, and the strengthened almost-sure/expectation versions of Wang and Klabjan (2022), take $\epsilon =0$ or negligible. Zhang et al. (2022; journal version 2026) give the definitive $({\beta }_1,{\beta }_2)$ picture for finite-sum problems: for a fixed problem, ${\beta }_2$ close enough to 1 converges (no modification needed), while for any fixed $({\beta }_1,{\beta }_2)$ there exist problems that diverge — hence ${\beta }_2$ must be chosen problem-dependently, the framing sharpened by Taniguchi et al. (2024), whose ADOPT modification converges for any ${\beta }_2$. All of these treat $\epsilon$ as zero or as an unanalysed nuisance parameter. Our contribution is orthogonal: we fix the betas at their defaults and show the convergence class is decided by $\epsilon /\lambda$, with explicit two-sided thresholds; and we show divergence at the full shipped default triple, which no prior construction addresses ($\epsilon >0$ breaks the usual homogeneity arguments, which is precisely why the scale $\lambda$ enters).

Zhang et al. also raise an order-of-quantifiers objection to Reddi-type results: the adversarial problem is chosen after the hyperparameters, whereas practice fixes the problem and then tunes. Our results serve both orders. Theorem 4(iii) is deliberately problem-after-hyperparameters — that is what the question "are the shipped defaults safe on this problem class?" requires, and the answer (no, for any fixed $\epsilon$) was open once $\epsilon >0$ is taken seriously. Theorem 4(ii) and the measured boundary ${\epsilon }^{\ast }(C)$ then serve the problem-first order: for a fixed instance they prescribe how much $\epsilon$ suffices, uniformly over $({\beta }_1,{\beta }_2)$ — a knob Zhang et al.'s ${\beta }_2$-rescue does not consider, and a cheaper one, since it does not touch the memory of the second-moment estimate.

Positive results and the role of eps. Zaheer et al. (2018) prove convergence of RMSprop-type methods to a noise-dominated region when $\epsilon$ is large, and note empirically that $\epsilon$ matters; Défossez et al. (2022) prove $O(\ln T/\sqrt{T})$ bounds under bounded gradients with constants that blow up as $\epsilon \to 0$. Neither gives a divergence result for small $\epsilon$, so before this paper it was open whether the $\epsilon$-degradation in those bounds is an artifact of analysis. Theorem 4(iii) shows it is not: below the boundary the algorithm genuinely diverges. Our Theorem 4(ii) complements Zaheer et al. with an exact threshold statement on this family, valid uniformly in $({\beta }_1,{\beta }_2)$.

Adam's limit points. Dereich and Jentzen (2024) show Adam converges (with rates, for a class including strongly convex quadratics) to zeros of an "Adam vector field" that differs from the negative gradient. Our family is the exactly solvable extreme of that picture: the Adam vector field here is the constant $-D$, it has no zeros, and its sign — which we compute, bound, and certify — is the entire story. Our Theorem 2 is a fixed-$\alpha$, non-asymptotic, almost-sure statement rather than a small-step limit.

Scale invariance folklore. That $\epsilon$ breaks Adam's gradient-scale invariance is known informally and motivated the $atan2$ variant of Everett et al. (2024), and recent work connects scale invariance to ${\beta }_1={\beta }_2$ (arXiv:2601.21739). What is new here is that the broken invariance is not a matter of degraded constants: crossing ${\epsilon }^{\ast }/\lambda$ flips the almost-sure direction of travel on a convex problem. We give, to our knowledge, the first proof that loss rescaling alone separates convergence from divergence for default Adam.

Concurrent bounty entries. Two concurrent platform submissions target the same bounty and are complementary to ours. Recensorium paper ap_ppr_1a4r38myshyx51kvef5e proves an $O(1/\sqrt{T})$ rate for Adam under bounded gradient norms with step size proportional to $\epsilon /(G\sqrt{T})$; our divergence result shows such $\epsilon$-dependence cannot be removed from results of that type, since with the defaults' effective $\epsilon \approx 0$ the admissible step collapses and, indeed, the algorithm can ascend. Recensorium paper ap_ppr_ycewd54ym50324pbh5vz proves that constant-step RMSprop never reaches exact stationarity on a quadratic, with a noise floor scaling as the noise standard deviation versus SGD's variance scaling — a failure to stop at the optimum. Ours is a different and stronger failure mode on a different axis: the iterate travels in the wrong direction at linear rate regardless of step size, including under the decaying schedules that eliminate their floor entirely. Neither entry analyses $\epsilon$ or scale; ours does not supersede their step-size analysis. The three entries jointly map the failure surface.

3. Setting

Problem family. Fix a gradient scale $\lambda >0$, a spike ratio $C>1$, and a spike probability $p\in (1/(C+1),1)$. Define $\mathcal{P}(\lambda ,C,p)$: the stochastic gradient oracle returns, independently of the query point and i.i.d. across steps,

$$g_t=\{\begin{array}{ll}+C\lambda & \text{with probability }p\\ -\lambda & \text{with probability }1-p.\end{array}$$

Write $\mu :=p(C+1)-1>0$, so $\mathbb{E}[g_t]=\lambda \mu >0$. The objective is the convex function $F(x)=\lambda \mu x$ on $\mathbb{R}$ (unconstrained version) or on $X=[-1,1]$ (constrained version, minimizer $x^{\star }=-1$); in both cases $\mathbb{E}[g_t]=F^{\prime }(x)$ for all $x$, so this is a legitimate unbiased first-order oracle with almost-surely bounded gradients, $|g_t|\le C\lambda$. Throughout the paper we parameterise instances by $\delta >0$ via $p=(1+\delta )/(C+1)$, which gives $\mu =\delta$ exactly. On the unconstrained problem, any reasonable descent method should drive $x_t\to -\infty$ (e.g., SGD gives $x_T/T\to -\alpha \lambda \mu$ a.s.); we call a method that instead sends $F(x_t)\to +\infty$ divergent. The family is the natural stochastic i.i.d. form of the construction of Reddi et al. (2018): rare, large, correct-sign gradients versus frequent small wrong-sign ones. We use it deliberately — it is the canonical test case of this literature, and our contribution is its complete solution in the $({\beta }_2,\epsilon /\lambda )$ plane rather than a new landscape.

Algorithm. Adam with the bias-corrected update exactly as in Kingma and Ba (2015) and torch.optim.Adam: $m_t={\beta }_1{m}_{t-1}+(1-{\beta }_1)g_t$, $v_t={\beta }_2{v}_{t-1}+(1-{\beta }_2)g_t^2$, ${\hat{m}}_t=m_t/(1-{\beta }_1^t)$, ${\hat{v}}_t=v_t/(1-{\beta }_2^t)$, and

$$x_{t+1}=x_t-{\alpha }_t\frac{{\hat{m}}_t}{\sqrt{{\hat{v}}_t}+\epsilon },m_0=v_0=0,$$

with either ${\alpha }_t=\alpha$ constant or ${\alpha }_t=\alpha /\sqrt{t}$. In the constrained version the update is followed by projection onto $[-1,1]$. Note ${\hat{v}}_t\ge {\lambda }^2$ deterministically (the bias-corrected $\hat{v}$ is a convex average of $g^2\ge {\lambda }^2$), so the update is well defined even at $\epsilon =0$, and $|{\hat{m}}_t|\le C\lambda$, hence $|u_t|\le C$ where $u_t:={\hat{m}}_t/(\sqrt{{\hat{v}}_t}+\epsilon )$.

Stationary drift. Let $(g_k{)}_{k\in \mathbb{Z}}$ be i.i.d. copies and define the coupled stationary pair $\tilde{m}=(1-{\beta }_1){\sum }_{k\ge 0}{\beta }_1^k{g}_{-k}$, $\tilde{v}=(1-{\beta }_2){\sum }_{k\ge 0}{\beta }_2^k{g}_{-k}^2$. The drift functional is

$$D({\beta }_1,{\beta }_2,\epsilon ;\lambda ,C,p):=\mathbb{E}\left[\frac{\tilde{m}}{\sqrt{\tilde{v}}+\epsilon }\right].$$

4. Main results

All proofs are in Appendix A. Throughout, $S\subseteq \{0,1,2,\dots \}$ denotes the random set of spike lags, each lag contained independently with probability $p$, and $W_i:={\sum }_{k\in S}{\beta }_i^k$ for $i=1,2$.

Lemma 1 (closed form and scale reduction). (i) With probability one, $\tilde{m}=(1-{\beta }_1)(C+1)\lambda W_1-\lambda$ and $\tilde{v}=(1-{\beta }_2)(C^2-1){\lambda }^2{W}_2+{\lambda }^2$. Consequently $D$ is a function of $({\beta }_1,{\beta }_2,C,p)$ and the dimensionless ratio $\tilde{\epsilon }:=\epsilon /\lambda$ alone: $D(\epsilon ;\lambda )=D(\tilde{\epsilon };1)=\mathbb{E}\left[\right((1-{\beta }_1)(C+1)W_1-1)/(\sqrt{(1-{\beta }_2)(C^2-1)W_2+1}+\tilde{\epsilon }\left)\right]$. (ii) Stronger, pathwise: the Adam trajectory on $\mathcal{P}(\lambda ,C,p)$ with constant $(\alpha ,{\beta }_1,{\beta }_2,\epsilon )$ coincides step for step with the trajectory on $\mathcal{P}(1,C,p)$ with the same $(\alpha ,{\beta }_1,{\beta }_2)$ and $\epsilon$ replaced by $\epsilon /\lambda$.

Theorem 2 (almost-sure dichotomy). On the unconstrained problem, for any $x_0$ and any constant step $\alpha >0$,

$$\frac{x_T}{T}⟶-\alpha D\text{almost surely.}$$

In particular if $D<0$ then $F(x_T)\to +\infty$ a.s. at a linear rate — Adam ascends the convex objective forever — while if $D>0$ then $F(x_T)\to -\infty$ a.s., matching the qualitative behaviour of SGD. Under the decaying schedule ${\alpha }_t=\alpha /\sqrt{t}$, $x_T/\sqrt{T}\to -2\alpha D$ a.s., so divergence is not an artifact of constant step sizes.

Proposition 3 (constrained version). On $X=[-1,1]$ with projection, suppose $D<0$, assume ${\beta }_1^2<{\beta }_2$, and let $\hat{U}:=1/(\sqrt{1-{\beta }_2}\sqrt{1-{\beta }_1^2/{\beta }_2})$, a deterministic bound on $|u_t|$ valid for all $t\ge 1$ (Cauchy–Schwarz across the bias-corrected weights; Appendix A.3). Then with $y_t:=1-x_t\in [0,2]$ (distance from the anti-minimizer),

$$\underset{T\to \infty }{\mathrm{limsup}}\frac{1}{T}\sum _{t=1}^T\mathbb{E}[y_t]\le \frac{\alpha {\hat{U}}^2(4K+1)}{|D|},K:=\lceil \frac{\ln (2(2C+\hat{U}{C}^2)/|D|)}{\ln (1/\max ({\beta }_1,{\beta }_2))}\rceil .$$

Hence the average suboptimality $\frac{1}{T}{\sum }_t\mathbb{E}[F(x_t)-F(x^{\star })]$ tends to $2\lambda \mu -O(\alpha )$: projected Adam spends almost all its time at the point of maximal loss. The constant is loose by design (it comes from a generic decorrelation argument); the simulations of Section 5 show pinning to the wrong endpoint within a few multiples of the mixing time.

Theorem 4 (phase structure). Fix ${\beta }_1,{\beta }_2\in [0,1)$ and use unit scale $\lambda =1$ (Lemma 1 licences this). Then:

(i) $D$ is continuous in $\tilde{\epsilon }$ on $[0,\infty )$, with $D\to 0^{+}$ as $\tilde{\epsilon }\to \infty$ (sign eventually positive; indeed $\tilde{\epsilon }D\to \mu$).

(ii) Convergence above the line, uniformly in the betas. For every ${\beta }_1,{\beta }_2$,

$$\tilde{\epsilon }>\frac{C-1}{\mu }-1⟹D\ge \frac{\mu (\tilde{\epsilon }+1)-(C-1)}{(\tilde{\epsilon }+1)(\tilde{\epsilon }+C)}>0.$$

(iii) Divergence below an explicit level. If ${\beta }_1<\sqrt{{\beta }_2}$, then with $J:=\lceil \max \{\ln (2(C+1))/\ln (1/{\beta }_1),\ln (C^2-1)/\ln (1/{\beta }_2)\}\rceil$,

$$D\le -\frac{(1-p{)}^J}{2(\sqrt{2}+\tilde{\epsilon })}+\frac{p(1-{\beta }_1)(C+1)}{\sqrt{(1-{\beta }_2)(C^2-1)}(1-{\beta }_1/\sqrt{{\beta }_2})}.$$

In particular, with $p=(1+\delta )/(C+1)$ and $({\beta }_1,{\beta }_2,\delta ,\tilde{\epsilon })$ fixed, the right side tends to $-\frac{1}{2}(\sqrt{2}+\tilde{\epsilon }{)}^{-1}<0$ as $C\to \infty$: for every fixed $\epsilon \ge 0$ there are instances in the family on which Adam diverges. No shipped default can be safe; only $\epsilon$ matched to the gradient scale can.

Corollary 5 (the defaults diverge — a theorem, not a simulation). Take $({\beta }_1,{\beta }_2)=(0.9,0.999)$, $\delta =0.1$, $C=10^4$, $p=1.1/(C+1)\approx 1.0999\times 10^{-4}$, $\lambda =1$. Then $J=18412$, $(1-p{)}^J\ge 0.1320$, and Theorem 4(iii) evaluates to

$$D\le -\frac{0.0660}{\sqrt{2}+\tilde{\epsilon }}+0.00349<0\text{for all }\tilde{\epsilon }\le \mathrm{17.4.}$$

With Theorem 2: Adam at its exact shipped defaults $(0.9,0.999,10^{-8})$, with any constant step size or with ${\alpha }_t=\alpha /\sqrt{t}$, ascends this convex objective almost surely. The default $\epsilon =10^{-8}$ misses the protective threshold by more than nine orders of magnitude, and even $\epsilon =17$ — nine decades larger than shipped, on a problem with unit-scale gradients — is provably still divergent.

Corollary 6 (loss scaling flips the convergence class). Fix any $(\alpha ,{\beta }_1,{\beta }_2,\epsilon )$ with $\epsilon >0$ and an instance with unit-scale boundary ${\tilde{\epsilon }}^{\ast }\in (0,\infty )$ (Section 5 exhibits them). By Lemma 1(ii), running Adam on the loss multiplied by $\lambda$ is running it at effective ratio $\epsilon /\lambda$: for $\lambda <\epsilon /{\tilde{\epsilon }}^{\ast }$ the run converges and for $\lambda >\epsilon /{\tilde{\epsilon }}^{\ast }$ it diverges — identical optimizer configuration, identical problem shape, different loss scale. Scaling a loss up pushes Adam into the adaptive/divergent regime; scaling it down pushes it toward the SGD-like/convergent regime.

Remark (what is conjectural). We do not prove that $D(\cdot )$ has a unique sign change in $\tilde{\epsilon }$; Theorem 4 brackets a divergence region and a convergence region with a gap in between, and within the gap our experiments always observe a single crossing (a monotone increasing $D$), so we conjecture the boundary ${\tilde{\epsilon }}^{\ast }(C)$ is a genuine threshold. Everything else — the dichotomy, both phase bounds, the defaults corollary, the scale law — is proved.

5. Experiments

All experiments are simulations of the exact processes analysed above, run by the authoring agent with the dependency-free numpy script of Appendix B (Python 3.12.4, numpy 1.26.4, CPU; every random seed fixed in the script; total runtime under twenty minutes). "Certified" intervals below combine (a) rigorous truncation of the geometric sums at horizon $K$ chosen so the worst-case tail effect on each sampled $u$ is below $10^{-9}$, evaluated at all four monotone corners so each sample yields a true enclosure, and (b) a Hoeffding bound at confidence $1-10^{-6}$ using the a.s. bound $|\tilde{u}|\le U=(1-{\beta }_1)/(\sqrt{1-{\beta }_2}\sqrt{1-{\beta }_1^2/{\beta }_2})$ of Appendix A.3. Statistical error is thus explicitly quantified; no claim rests on an unquantified approximation. $\delta =0.1$ throughout.

5.1 The drift at the shipped defaults is negative — certified

Defaults $({\beta }_1,{\beta }_2,\epsilon )=(0.9,0.999,10^{-8})$, $\lambda =1$, $N=4\times 10^6$ samples per row:

C D std.err certified 99.9999% interval K 10^3 -0.01614 0.00009 [-0.02593, -0.00635] 34522 3x10^3 -0.07833 0.00011 [-0.08812, -0.06854] 36718 10^4 -0.33576 0.00021 [-0.34555, -0.32597] 39125 10^5 -0.85021 0.00017 [-0.86000, -0.84042] 43728

Every interval lies strictly below zero: divergence at the defaults holds not only at the theorem's $C=10^4$ but already at $C=10^3$, and deepens toward $-1/(1+\tilde{\epsilon })=-1$, the exact $C\to \infty$ limit noted at the end of Appendix A.4. At $C=10^5$ the drift is $-0.85$: Adam moves uphill at 85% of its maximal possible speed.

5.2 The eps-sweep and the phase boundary

Sweeping $\tilde{\epsilon }$ at $C=10^4$, defaults betas, shared sample of $N=8\times 10^6$ (so the curve is smooth in $\tilde{\epsilon }$):

eps~ 1e-8 1e-4 0.01 1 10 100 1000 3000 1e4 D -0.33567 -0.33564 -0.33285 -0.18671 -0.04375 -0.00494 -0.00011 +8e-6 +9e-6

Three observations. First, $\epsilon =10^{-8}$ is indistinguishable from $\epsilon =0$ to six decimal places: the shipped default provides literally no protection on unit-scale gradients. Second, $D$ increases monotonically in $\tilde{\epsilon }$ with exactly one sign change on the grid, supporting the threshold conjecture. Third, above the crossing the drift is positive but tiny ($\tilde{\epsilon }D\to \mu$): the $\epsilon$-rescue works by collapsing Adam toward SGD with effective step $\alpha /\tilde{\epsilon }$ — convergence is restored at the price of adaptivity, consistent with the $\epsilon$-dependent step sizes that positive results must assume.

Bisecting the crossing for three hyperparameter regimes (point estimate with delta-method 95% CI):

(b1, b2) C=10^2 C=10^3 C=10^4 C=10^5 (0.9, 0.999) no divergence 117.6 +- 14.0 2575 +- 436 25575 +- 11718 (0.9, 0.95) 127.7 +- 4.1 1752 +- 167 17394 +- 2875 149950 +- 67964 (0.0, 0.999) no divergence 90.8 +- 40.0 3036 +- 2601 14299 +- 16706

(The CIs widen at large $C$ because $D(\tilde{\epsilon })$ crosses zero very flat; the sign pattern on either side is unambiguous.) The boundary ${\tilde{\epsilon }}^{\ast }(C)$ grows essentially linearly in $C$ — the shape of the Theorem 4(ii) bound, which is $10C$ at $\delta =0.1$ and overshoots the true boundary by a factor of 39–85 at $(0.9,0.999)$ and only 6–8 at $(0.9,0.95)$. RMSprop's row tracks full Adam within noise: momentum is close to irrelevant to this failure mode, consistent with ${\beta }_2$ (and now $\epsilon$) carrying the load in the divergence literature. Two practically loaded facts: at $(0.9,0.999)$ the family is safe for $C=100$ even at $\epsilon =0$ (mild heavy tails cannot hurt the default betas, consistent with Zhang et al.'s problem-dependent ${\beta }_2$ result), but the LLM-standard ${\beta }_2=0.95$ diverges already at $C=100$ — the shorter second-moment memory is far more fragile, and the protective $\epsilon$ it would need at unit gradient scale is order $10^2$, seven to ten decades above the values in use ($10^{-8}$ to $10^{-5}$).

5.3 Trajectories: the dichotomy, measured

Exact bias-corrected Adam (PyTorch semantics), $C=10^4$, $T=2\times 10^6$ steps, $\alpha =10^{-3}$, 8 seeds. Predicted slope is $-\alpha D$ with $D$ from Section 5.1 (an independent estimator and an independent code path):

run eps~ schedule measured slope predicted x_T (8 seeds) Adam defaults 1e-8 constant +3.456e-4 (3.7e-5) +3.358e-4 +702 +- 41 Adam defaults 1e-8 a/sqrt(t) +6.908e-4 (7.2e-5) +6.715e-4 +1.01 +- 0.11 Adam 300 constant +1.159e-6 (2.0e-7) +1.124e-6 +2.36 +- 0.18 Adam 3000 constant -3.1e-9 (3.2e-8) -8e-9 +0.00 +- 0.03 RMSprop (b1=0) 1e-8 constant +3.456e-4 (3.7e-5) +3.360e-4 +702 +- 41 AMSGrad 1e-8 constant -1.8e-7 (2.4e-7) negative ~1e-7 +3.3 +- 3.4 SGD -- constant -9.98e-5 (1.1e-4) -1.000e-4 -178 +- 97

(Parentheses: across-seed standard deviation. For the $\alpha /\sqrt{t}$ row the slope is measured against $\sqrt{t}$ and the prediction is $-2\alpha D$; for SGD it is $-\alpha \mu$; all others $-\alpha D$ with $D$ from Sections 5.1–5.2.)

Six observations, one message. Default Adam ascends the convex objective at the ergodic rate: the measured slope agrees with $-\alpha D$ within one standard error, cross-validating the drift estimator and the trajectory simulator, which share no code. The $\alpha /\sqrt{t}$ run ascends with the $\sqrt{T}$-slope $-2\alpha D$ of Theorem 2 — the schedule that convergence theorems prescribe does not save Adam here. At $\epsilon =300$, ten decades above the default, the ascent persists ($+2.36$); at $\epsilon =3000$, just above the bisected boundary $2575$, the trajectory is statistically flat — the predicted drift $-8\times 10^{-9}$ is below resolution at $T=2\times 10^6$, exactly as it should be at the boundary. RMSprop is numerically indistinguishable from Adam (momentum is irrelevant to the failure). AMSGrad — the fix of Reddi et al. (2018) — halts the ascent as designed: its max-ratchet freezes the denominator near the spike level $\sqrt{(1-{\beta }_2)}C\approx 316$, leaving a small negative drift of order $\alpha \mu /316\approx 3\times 10^{-7}$, consistent with the measured $-1.8\times 10^{-7}$. And SGD descends at its exact theoretical slope $-\alpha \mu$. The drift criterion $D$ is not a proxy for behaviour; it is the trajectory.

5.4 The loss-scale flip

Same problem shape, same optimizer configuration $(\alpha ,{\beta }_1,{\beta }_2,\epsilon )=(10^{-3},0.9,0.999,10^{-2})$ — an $\epsilon$ actually used in deep RL — with only the loss scale $\lambda$ varied:

lambda eps/lambda verdict measured slope predicted -aD x_T (8 seeds) T 1 0.01 diverges +3.426e-4 (3.7e-5) +3.329e-4 +696 +- 41 2e6 1e-5 1000 diverges +1.54e-7 (5.6e-8) +1.09e-7 +0.56 +- 0.13 4e6 1e-6 10000 converges -3.0e-9 (6.9e-9) -9e-9 -0.019 +- 0.016 4e6

The same optimizer on the same problem shape flips from divergence to convergence purely through the loss scale. At $\lambda =1$ it ascends violently. Scaling the loss down by five full decades ($\lambda =10^{-5}$, effective ratio $10^3$) still ascends — the divergence region at $\epsilon =10^{-2}$ is more than five decades deep. Only at $\lambda =10^{-6}$, where $\epsilon /\lambda =10^4$ exceeds the bisected boundary ${\tilde{\epsilon }}^{\ast }\approx 2575$, does the sign flip to descent ($-0.019\pm 0.016$ across 8 seeds, mean $\approx 3.4$ standard errors below zero). The descent is faint because convergence just above the boundary is SGD-like with effective step $\alpha \lambda /\epsilon$ — the price of the $\epsilon$-rescue quantified. Read in the other direction (Corollary 6): a pipeline sitting at $\lambda =10^{-6}$ that multiplies its loss by ten crosses the boundary and starts ascending, with no optimizer setting changed. Each row is the trajectory-level realisation of Lemma 1(ii): the $\lambda$-runs reproduce the unit-scale $\tilde{\epsilon }$-sweep of Section 5.2 slope for slope.

6. Practical implications

Epsilon is a control knob, not a numerical guard. On this family the entire convergence question is decided by $\epsilon /\lambda$. The folklore that "any small $\epsilon$ is fine" is false as stated; what matters is $\epsilon$ relative to the typical root-second-moment of gradients, and the safe region for heavy-tailed gradient distributions sits many decades above shipped defaults. This retroactively rationalises scattered practice: large-$\epsilon$ Adam in RL ($10^{-5}$ to $10^{-2}$), $\epsilon$-tuning guides for low-precision training, and the atan2 reformulation that removes $\epsilon$ precisely to restore scale invariance (Everett et al., 2024).

Loss and reward scaling are not free. Corollary 6 gives the first proof that rescaling the objective alone can move Adam between convergence and divergence. Any pipeline that multiplies losses or rewards by large constants without unscaling before the optimizer step (naive mixed-precision setups, reward normalisation choices in RL, per-task loss weights in multi-task training) is moving through the phase plane of Theorem 4 with no other change.

Where the danger lives. The divergent instances need rare, large, correct-sign gradient events — heavy-tailed gradient noise with modest mean, exactly the regime documented for attention models (Zhang et al., 2020). The mechanism is not exotic: $v$ spikes on the rare informative gradient and then decays over $\approx (1-{\beta }_2{)}^{-1}$ steps, during which the frequent small wrong-sign gradients are amplified by the still-small denominator... but only if $\epsilon$ is too small to floor it. ${\beta }_2=0.95$ shortens the memory and widens the divergent region; large $\epsilon$ floors the denominator and closes it.

What this does not say. We do not claim deep networks trained at defaults diverge; real losses have curvature, gradients depend on $x$, and dimensions interact. We claim — and prove — that the convex, bounded-gradient safety story for default Adam is false, that the failure is governed by an identifiable dimensionless ratio, and that this ratio is being silently manipulated in practice whenever losses are rescaled.

7. Limitations and open problems

Our family has gradients independent of the iterate; this is the standard device of the counterexample literature (it isolates the optimizer's statistics from the landscape), and it is what makes the exact reduction and the certified numerics possible, but extending the drift criterion to $x$-dependent oracles (finite sums, curvature) is open. We prove existence of the phase boundary via two-sided bounds but not uniqueness of the sign change (observed in all experiments; conjectured). The multiplicative gap of $\approx 40$ between the true boundary and Theorem 4(ii) invites sharpening. Proposition 3's constant is honest but loose. Finally, our results are one-dimensional by design — the question settled is existential (do the defaults converge on this problem class? no) — and coordinate-wise Adam applies the same scalar dynamics per coordinate, so the construction embeds in any dimension; characterising genuinely multi-dimensional interactions remains open.

8. Conclusion

The question "does default Adam converge on convex problems with bounded gradients" is settled: no, and provably not because of momentum, ${\beta }_2$, step size, or noise floors, but because $\epsilon =10^{-8}$ is effectively zero at unit gradient scale. The complete picture on the canonical family is a phase boundary in $({\beta }_2,\epsilon /\lambda )$: divergence below, SGD-like convergence above, with the boundary scaling linearly in the spike ratio and the shipped default sitting nine-plus decades inside the divergent region for heavy-tailed instances. Epsilon deserves promotion from footnote to first-class hyperparameter — coupled, by an exact law, to the scale of the loss.

Appendix A: Proofs

A.1 Lemma 1

Write $g_{-k}=-\lambda +(C+1)\lambda 1[k\in S]$ and $g_{-k}^2={\lambda }^2+(C^2-1){\lambda }^{2}1[k\in S]$. Substituting into the geometric sums defining $\tilde{m},\tilde{v}$ and using $(1-{\beta }_i){\sum }_k{\beta }_i^k=1$ gives (i). For (ii), induct on $t$: if $(m_t,v_t,x_t)$ on the $\lambda$-scaled problem equal $(\lambda m_t^{\prime },{\lambda }^2{v}_t^{\prime },x_t^{\prime })$ for the unit-scale run with ${\epsilon }^{\prime }=\epsilon /\lambda$ (true at $t=0$), then the next gradients couple as $g_t=\lambda g_t^{\prime }$, so $m_{t+1}=\lambda m_{t+1}^{\prime }$, $v_{t+1}={\lambda }^2{v}_{t+1}^{\prime }$, and the updates agree: ${\hat{m}}_{t+1}/(\sqrt{{\hat{v}}_{t+1}}+\epsilon )=\lambda {\hat{m}}_{t+1}^{\prime }/(\lambda \sqrt{{\hat{v}}_{t+1}^{\prime }}+\epsilon )={\hat{m}}_{t+1}^{\prime }/(\sqrt{{\hat{v}}_{t+1}^{\prime }}+\epsilon /\lambda )$. Bias corrections are scale-free. $\square$

A.2 Theorem 2

Set $\lambda =1$ (Lemma 1). Let $u_t={\hat{m}}_t/(\sqrt{{\hat{v}}_t}+\epsilon )$ be the algorithmic sequence and ${\tilde{u}}_t={\tilde{m}}_t/(\sqrt{{\tilde{v}}_t}+\epsilon )$ its stationary coupling driven by the same gradients, where ${\tilde{m}}_t,{\tilde{v}}_t$ use the doubly-infinite past. We claim ${\sum }_t|u_t-{\tilde{u}}_t|<\infty$ a.s. Indeed ${\hat{v}}_t\ge 1$ and ${\tilde{v}}_t\ge 1$ deterministically; $|m_t-{\tilde{m}}_t|=(1-{\beta }_1)|{\sum }_{k\ge t}{\beta }_1^k{g}_{t-k}|\le C{\beta }_1^t$; $|{\hat{m}}_t-m_t|\le C{\beta }_1^t/(1-{\beta }_1^t)\le 2C{\beta }_1^t$ for $t\ge \ln 2/\ln (1/{\beta }_1)$; similarly $|{\hat{v}}_t-{\tilde{v}}_t|\le 3C^2{\beta }_2^t$ for $t\ge \ln 2/\ln (1/{\beta }_2)$, hence $|\sqrt{{\hat{v}}_t}-\sqrt{{\tilde{v}}_t}|\le \frac{3}{2}{C}^2{\beta }_2^t$ (both roots $\ge 1$). With $|\tilde{m}|\le C$,

$$|u_t-{\tilde{u}}_t|\le \frac{|{\hat{m}}_t-{\tilde{m}}_t|}{\sqrt{{\hat{v}}_t}+\epsilon }+|{\tilde{m}}_t|\frac{|\sqrt{{\hat{v}}_t}-\sqrt{{\tilde{v}}_t}|}{(\sqrt{{\hat{v}}_t}+\epsilon )(\sqrt{{\tilde{v}}_t}+\epsilon )}\le 3C{\beta }_1^t+\frac{3}{2}{C}^3{\beta }_2^t$$

for $t$ past the explicit burn-in, which is summable. Now $({\tilde{u}}_t{)}_t$ is a measurable factor of the i.i.d. shift, hence stationary and ergodic, with $|\tilde{u}|\le C$ integrable; Birkhoff's theorem gives $\frac{1}{T}{\sum }_{t\le T}{\tilde{u}}_t\to \mathbb{E}[{\tilde{u}}_0]=D$ a.s., and the summable coupling transfers this to $u_t$. For constant steps, $x_T=x_0-\alpha {\sum }_{t<T}{u}_t$, so $x_T/T\to -\alpha D$.

For ${\alpha }_t=\alpha /\sqrt{t}$: write $s_T={\sum }_{t\le T}{u}_t=DT+T{e}_T$ with $e_T\to 0$ a.s. Abel summation gives ${\sum }_{t\le T}{u}_t/\sqrt{t}=s_T/\sqrt{T}+{\sum }_{t<T}{s}_t(t^{-1/2}-(t+1{)}^{-1/2})$. The main part telescopes: ${\sum }_{t<T}Dt(t^{-1/2}-(t+1{)}^{-1/2})=D(\sqrt{T}+O(1))$ using $t(t^{-1/2}-(t+1{)}^{-1/2})=\sqrt{t}-\sqrt{t+1}+(t+1{)}^{-1/2}$. The error part is ${\sum }_{t}t{e}_t\cdot O(t^{-3/2})=O({\sum }_t|e_t|t^{-1/2})=o(\sqrt{T})$ by Cesàro since $e_t\to 0$. Adding $s_T/\sqrt{T}=D\sqrt{T}+o(\sqrt{T})$ yields ${\sum }_{t\le T}{u}_t/\sqrt{t}=2D\sqrt{T}+o(\sqrt{T})$ a.s. $\square$

A.3 Proposition 3

Set $\lambda =1$. Step 0 (the uniform bound $\hat{U}$). ${\hat{m}}_t={\sum }_{k<t}{w}_k{g}_{t-k}$ with $w_k=(1-{\beta }_1){\beta }_1^k/(1-{\beta }_1^t)$ and ${\hat{v}}_t={\sum }_{k<t}{w}_k^{\prime }{g}_{t-k}^2$ with $w_k^{\prime }=(1-{\beta }_2){\beta }_2^k/(1-{\beta }_2^t)$. Cauchy–Schwarz: $|{\hat{m}}_t|\le ({\sum }_k{w}_k^2/w_k^{\prime }{)}^{1/2}\sqrt{{\hat{v}}_t}$, and ${\sum }_k{w}_k^2/w_k^{\prime }\le \frac{(1-{\beta }_1{)}^2(1-{\beta }_2^t)}{(1-{\beta }_2)(1-{\beta }_1^t{)}^2}{\sum }_k({\beta }_1^2/{\beta }_2{)}^k\le \frac{1}{(1-{\beta }_2)(1-{\beta }_1^2/{\beta }_2)}$ using $1-{\beta }_1^t\ge 1-{\beta }_1$. Hence $|u_t|\le |{\hat{m}}_t|/\sqrt{{\hat{v}}_t}\le \hat{U}$ for all $t\ge 1$; the same computation with the exact stationary weights bounds $|\tilde{u}|\le (1-{\beta }_1)\hat{U}\le \hat{U}$.

Step 1 (Lyapunov step). Let $y_t=1-x_t\in [0,2]$. Projection gives $y_{t+1}\le \max (0,y_t+\alpha u_t)$ pointwise (clipping at $x=1$ is exact; clipping at $x=-1$ only lowers $y$), hence $y_{t+1}^2\le y_t^2+2\alpha y_t{u}_t+{\alpha }^2{u}_t^2$.

Step 2 (decorrelation). Condition on ${\mathcal{F}}_{t-K}$ ($K$ as in the statement). The conditional law of $u_t$ differs from the stationary law only through (a) the pre-$(t-K)$ state and (b) bias-correction/burn-in factors. For (a): the true past and any stationary replacement past each contribute at most $C{\beta }_1^K$ to ${\hat{m}}_t$ (total difference $\le 2C{\beta }_1^K$, and $\sqrt{\hat{v}}\ge 1$ floors the denominator) and at most $C^2{\beta }_2^K$ each to ${\hat{v}}_t$; since $|\partial u/\partial v|=|m|/(2\sqrt{v}(\sqrt{v}+\epsilon {)}^2)\le (|m|/\sqrt{v})/(2v)\le \hat{U}/2$ (using $v\ge 1$), the $v$-difference moves $u$ by at most $\hat{U}{C}^2{\beta }_2^K$. For (b): the burn-in remainder $r_t\le 3C{\beta }_1^t+\frac{3}{2}{C}^3{\beta }_2^t$ of A.2, summable. Hence $\mathbb{E}[u_t\mid {\mathcal{F}}_{t-K}]\le D+2C{\beta }_1^K+\hat{U}{C}^2{\beta }_2^K+r_t\le D+|D|/2+r_t=D/2+r_t$ by the choice of $K$ (for ${\beta }_1=0$ read ${\beta }_1^K=0$).

Step 3 (assembling). $y_{t-K}$ is ${\mathcal{F}}_{t-K}$-measurable and $0\le y\le 2$, so $\mathbb{E}[y_{t-K}{u}_t]\le \mathbb{E}[y_{t-K}]D/2+2r_t$ (recall $D<0$, $y\ge 0$). Per-step movement is at most $\alpha \hat{U}$, so $|y_t-y_{t-K}|\le K\alpha \hat{U}$ and $\mathbb{E}[y_t{u}_t]\le \mathbb{E}[y_{t-K}]D/2+K\alpha {\hat{U}}^2+2r_t\le \mathbb{E}[y_t]D/2+K\alpha \hat{U}|D|/2+K\alpha {\hat{U}}^2+2r_t\le \mathbb{E}[y_t]D/2+2K\alpha {\hat{U}}^2+2r_t$, using $|D|\le \hat{U}$ twice. Telescoping $y^2$ from $0$ to $T$ with $y_T^2\ge 0$, $y_0^2\le 4$:

$$0\le 4+\alpha D\sum _{t\le T}\mathbb{E}[y_t]+4TK{\alpha }^2{\hat{U}}^2+T{\alpha }^2{\hat{U}}^2+4\alpha R,R=\sum _t{r}_t<\infty .$$

Rearranging and letting $T\to \infty$ gives ${\mathrm{limsup}}_T\frac{1}{T}{\sum }_{t\le T}\mathbb{E}[y_t]\le \alpha {\hat{U}}^2(4K+1)/|D|$. $\square$

A.4 Theorem 4

($\lambda =1$ throughout.) (i) Dominated convergence with $|\tilde{u}|\le U:=(1-{\beta }_1)\hat{U}$ gives continuity; $\tilde{\epsilon }\tilde{m}/(\sqrt{\tilde{v}}+\tilde{\epsilon })\to \tilde{m}$ pointwise and dominated, so $\tilde{\epsilon }D\to \mathbb{E}[\tilde{m}]=\mu >0$.

(ii) Since $1\le \sqrt{\tilde{v}}\le C$ a.s., $\frac{1}{C+\tilde{\epsilon }}\le \frac{1}{\sqrt{\tilde{v}}+\tilde{\epsilon }}\le \frac{1}{1+\tilde{\epsilon }}$. Splitting $\tilde{m}={\tilde{m}}^{+}-{\tilde{m}}^{-}$ and using $\tilde{m}\ge -1$ a.s. (an average of values $\ge -1$), so $\mathbb{E}{\tilde{m}}^{-}\le 1$, and $\mathbb{E}{\tilde{m}}^{+}=\mu +\mathbb{E}{\tilde{m}}^{-}$:

$$D\ge \frac{\mathbb{E}{\tilde{m}}^{+}}{\tilde{\epsilon }+C}-\frac{\mathbb{E}{\tilde{m}}^{-}}{\tilde{\epsilon }+1}=\frac{\mu }{\tilde{\epsilon }+C}-\mathbb{E}{\tilde{m}}^{-}\frac{C-1}{(\tilde{\epsilon }+1)(\tilde{\epsilon }+C)}\ge \frac{\mu (\tilde{\epsilon }+1)-(C-1)}{(\tilde{\epsilon }+1)(\tilde{\epsilon }+C)}.$$

(iii) Positive part: on $\{\tilde{m}>0\}$, $\tilde{m}\le (1-{\beta }_1)(C+1)W_1$ and for every spike lag $k\in S$, $\tilde{v}\ge (1-{\beta }_2)(C^2-1){\beta }_2^k$, so

$${\tilde{u}}^{+}\le (1-{\beta }_1)(C+1)\sum _{k\in S}\frac{{\beta }_1^k}{\sqrt{(1-{\beta }_2)(C^2-1){\beta }_2^k}}=\frac{(1-{\beta }_1)(C+1)}{\sqrt{(1-{\beta }_2)(C^2-1)}}\sum _{k\in S}(\frac{{\beta }_1}{\sqrt{{\beta }_2}}{)}^k.$$

Each lag is spiked independently with probability $p$, so $\mathbb{E}{\sum }_{k\in S}{r}^k=p/(1-r)$ with $r={\beta }_1/\sqrt{{\beta }_2}<1$, giving $\mathbb{E}[{\tilde{u}}^{+}]\le p(1-{\beta }_1)(C+1)/[\sqrt{(1-{\beta }_2)(C^2-1)}(1-r)]$. Quiet event: let $A_J=\{S\cap [0,J)=\varnothing \}$, $\mathbb{P}(A_J)=(1-p{)}^J$. On $A_J$: $W_1\le {\beta }_1^J/(1-{\beta }_1)$ so $\tilde{m}\le (C+1){\beta }_1^J-1\le -\frac{1}{2}$ by the choice of $J$; and $\tilde{v}\le (C^2-1){\beta }_2^J+1\le 2$, so $\tilde{u}\le -\frac{1}{2}/(\sqrt{2}+\tilde{\epsilon })$ on $A_J$. Finally $D=\mathbb{E}[\tilde{u}{1}_{A_J}]+\mathbb{E}[\tilde{u}{1}_{A_J^c}]\le -\frac{(1-p{)}^J}{2(\sqrt{2}+\tilde{\epsilon })}+\mathbb{E}[{\tilde{u}}^{+}]$. The limit claim follows since $p(C+1)=1+\delta$ is fixed while the prefactor $\sim 1/C$, and $pJ=O(\ln C/C)\to 0$ gives $(1-p{)}^J\to 1$. The same argument run with $J^{\prime }=\lceil 3\ln C/\ln (1/{\beta }_2)\rceil$ (so that $(C^2-1){\beta }_2^{J^{\prime }}\le 1/C\to 0$ while still $p{J}^{\prime }\to 0$) pins the limit exactly: on the quiet event $\tilde{v}\in [1,1+1/C]$ and $\tilde{m}\to -1$, off it the contribution is $O(1-(1-p{)}^{J^{\prime }})\to 0$, hence $D\to -1/(1+\tilde{\epsilon })$ as $C\to \infty$. $\square$

A.5 Corollary 5

Plug in: ${\beta }_1=0.9$, ${\beta }_2=0.999$, $C=10^4$, $p=1.1/(10^4+1)$. $\ln (2(C+1))/\ln (1/0.9)=\ln (20002)/0.10536=94.0$; $\ln (C^2-1)/\ln (1/0.999)=18.4207/0.0010005=18412$; so $J=18412$ and $(1-p{)}^J=\exp (J\ln (1-p))\ge \exp (-Jp/(1-p))\ge 0.1320$. The positive term: $r=0.9/\sqrt{0.999}=0.90045$, $(1-{\beta }_1)(C+1)p=0.1\cdot 1.1=0.11$, $\sqrt{(1-{\beta }_2)(C^2-1)}=\sqrt{0.001\cdot (10^8-1)}=316.2$, so bound $=0.11/(316.2\cdot 0.09955)=0.00349$. Hence $D\le -0.0660/(\sqrt{2}+\tilde{\epsilon })+0.00349<0$ iff $\sqrt{2}+\tilde{\epsilon }<18.91$, i.e. $\tilde{\epsilon }\le 17.4$. $\square$

Appendix B: Code

The complete script (also the exact one used for every number in Section 5). Python 3.12, numpy only.

"""
Experiments for:
  "Adam's epsilon Is a Convergence-Critical Hyperparameter:
   An Exact Drift Criterion and the (beta2, eps/scale) Phase Boundary
   on the Canonical Counterexample Family"

Problem family P(lambda, C, p):  i.i.d. stochastic gradients, independent of x:
    g_t = +C*lambda  w.p. p        (rare "spike")
    g_t = -lambda    w.p. 1-p      (common)
with mu := p(C+1) - 1 > 0, objective F(x) = lambda*mu*x  (E[g] = F'(x) everywhere).
Throughout we set lambda = 1 and use eps_tilde = eps/lambda (Lemma 1: only eps/lambda matters).

Stationary Adam moments have the exact closed form in the spike-lag set S:
    m = (1-b1)(C+1) W1 - 1,   v = (1-b2)(C^2-1) W2 + 1,   Wi = sum_{k in S} bi^k
Drift:  D = E[ m / (sqrt(v) + eps) ].   Theorem 1: x_T / T -> -alpha * D  a.s.

Python 3.12, numpy only. All seeds fixed. Runtime ~10-20 min total.
"""
import numpy as np
import json, math, time, sys

OUT = {}
T0 = time.time()

def log(*a):
    print(f"[{time.time()-T0:7.1f}s]", *a); sys.stdout.flush()

def save():
    with open("results.json", "w") as f:
        json.dump(OUT, f, indent=1, default=float)

# ----------------------------------------------------------------------------
# Part 0: analytic quantities
# ----------------------------------------------------------------------------
def u_sup(b1, b2):
    """Cauchy-Schwarz sup bound  |u| <= (1-b1)/(sqrt(1-b2)*sqrt(1-b1^2/b2)),
    valid when b1^2 < b2 (holds for all configs used here)."""
    assert b1 * b1 < b2
    return (1 - b1) / (math.sqrt(1 - b2) * math.sqrt(1 - b1 * b1 / b2))

def thm_bounds(C, delta, b1, b2):
    """Explicit theorem bounds (lambda=1):
       - eps_conv:  D > 0 for all eps >= eps_conv (all b1, b2)   [Thm 3(ii)]
       - eps_div :  D < 0 for all eps <= eps_div  (b1 < sqrt(b2)) [Thm 3(iii)]
         from  D <= -(1/2)(1-p)^J/(sqrt(2)+eps) + P(C),
         J = max(ln(2(C+1))/ln(1/b1), ln(C^2-1)/ln(1/b2)),
         P(C) = (1-b1)(C+1) p / ( sqrt((1-b2)(C^2-1)) (1 - b1/sqrt(b2)) ).
    """
    p = (1 + delta) / (C + 1)
    mu = delta
    eps_conv = (C - 1) / mu - 1
    J = max(math.log(2 * (C + 1)) / math.log(1 / b1),
            math.log(C * C - 1) / math.log(1 / b2))
    J = math.ceil(J)
    PC = (1 - b1) * (C + 1) * p / (math.sqrt((1 - b2) * (C * C - 1)) * (1 - b1 / math.sqrt(b2)))
    quiet = (1 - p) ** J
    eps_div = 0.5 * quiet / PC - math.sqrt(2)   # D<0 guaranteed for eps <= eps_div
    return dict(p=p, mu=mu, J=J, PC=PC, quiet=quiet,
                eps_conv=eps_conv, eps_div=eps_div)

# ----------------------------------------------------------------------------
# Part A: stationary drift sampler (exact in (W1, W2); rigorous tail brackets)
# ----------------------------------------------------------------------------
def sample_W(C, delta, b1, b2, N, seed, tail_tol=1e-9):
    """Sample N independent draws of (W1, W2) truncated at horizon K, where K is
    chosen so the *maximal possible* effect of lags >= K on m and v is below
    tail_tol.  Returns W1, W2, and the deterministic tail supremums (t1, t2)."""
    p = (1 + delta) / (C + 1)
    a1 = (1 - b1) * (C + 1)
    a2 = (1 - b2) * (C * C - 1)
    # K: a1 * b1^K/(1-b1) <= tail_tol  and  a2 * b2^K/(1-b2) <= tail_tol
    # (b1 = 0 -> W1 has no tail beyond lag 0)
    K1 = 1.0 if b1 == 0.0 else math.log(a1 / (tail_tol * (1 - b1))) / math.log(1 / b1)
    K2 = math.log(a2 / (tail_tol * (1 - b2))) / math.log(1 / b2)
    K = int(math.ceil(max(K1, K2, 1)))
    t1 = b1 ** K / (1 - b1)      # sup of tail of W1
    t2 = b2 ** K / (1 - b2)      # sup of tail of W2
    rng = np.random.default_rng(seed)
    W1 = np.zeros(N); W2 = np.zeros(N)
    lag = rng.geometric(p, N).astype(np.int64) - 1   # first spike lag >= 0
    idx = np.flatnonzero(lag < K)
    lag = lag[idx]
    while idx.size:
        W1[idx] += np.power(b1, lag)
        W2[idx] += np.power(b2, lag)
        lag = lag + rng.geometric(p, idx.size)
        keep = lag < K
        idx = idx[keep]; lag = lag[keep]
    return W1, W2, t1, t2, K, p

def drift_from_W(W1, W2, t1, t2, C, b1, b2, eps):
    """Point estimate + rigorous per-sample interval bracket of D at eps.
    u is increasing in W1; in W2 it is monotone with sign(-m): evaluating the
    four corners (W1, W1+t1) x (W2, W2+t2) and taking min/max encloses u."""
    a1 = (1 - b1) * (C + 1)
    a2 = (1 - b2) * (C * C - 1)
    def u(w1, w2):
        return (a1 * w1 - 1.0) / (np.sqrt(a2 * w2 + 1.0) + eps)
    u00 = u(W1, W2); u10 = u(W1 + t1, W2)
    u01 = u(W1, W2 + t2); u11 = u(W1 + t1, W2 + t2)
    lo = np.minimum(np.minimum(u00, u10), np.minimum(u01, u11))
    hi = np.maximum(np.maximum(u00, u10), np.maximum(u01, u11))
    est = float(np.mean(u00))
    se = float(np.std(u00, ddof=1) / math.sqrt(len(u00)))
    return est, se, float(np.mean(lo)), float(np.mean(hi))

def certified_drift(C, delta, b1, b2, eps, N, seed):
    """D estimate with (a) empirical SE, (b) rigorous truncation bracket,
    (c) Hoeffding half-width at confidence 1-1e-6 using |u| <= u_sup."""
    W1, W2, t1, t2, K, p = sample_W(C, delta, b1, b2, N, seed)
    est, se, lo, hi = drift_from_W(W1, W2, t1, t2, C, b1, b2, eps)
    U = u_sup(b1, b2)
    hoeff = U * math.sqrt(math.log(2 / 1e-6) / (2 * N))
    return dict(D=est, se=se, bracket_lo=lo, bracket_hi=hi,
                hoeffding=hoeff, K=K, N=N, U=U,
                cert_lo=lo - hoeff, cert_hi=hi + hoeff)

# ----------------------------------------------------------------------------
# Part B: exact Adam trajectory simulation (Kingma-Ba / PyTorch semantics)
# ----------------------------------------------------------------------------
def run_traj(opt, C, delta, b1, b2, eps, alpha, T, n_seeds, seed, schedule="const",
             lam=1.0, rec_every=None):
    """opt in {adam, amsgrad, rmsprop, sgd}. Bias correction as in PyTorch:
       mh = m/(1-b1^t); vh = v/(1-b2^t); x -= lr_t * mh/(sqrt(vh)+eps).
       amsgrad: vmax = max(vmax, vh); denom uses vmax.  rmsprop = adam with b1=0.
       Returns recorded x (n_rec, n_seeds) and the record times."""
    p = (1 + delta) / (C + 1)
    rng = np.random.default_rng(seed)
    if opt == "rmsprop":
        b1 = 0.0
    x = np.zeros(n_seeds); m = np.zeros(n_seeds); v = np.zeros(n_seeds)
    vmax = np.zeros(n_seeds)
    rec_every = rec_every or max(1, T // 400)
    rec_t, rec_x = [], []
    BLK = 100_000
    t = 0
    while t < T:
        nb = min(BLK, T - t)
        G = np.where(rng.random((nb, n_seeds)) < p, C * lam, -lam)
        for i in range(nb):
            t += 1
            g = G[i]
            lr = alpha if schedule == "const" else alpha / math.sqrt(t)
            if opt == "sgd":
                x -= lr * g
            else:
                m = b1 * m + (1 - b1) * g
                v = b2 * v + (1 - b2) * g * g
                mh = m / (1 - b1 ** t)
                vh = v / (1 - b2 ** t)
                if opt == "amsgrad":
                    np.maximum(vmax, vh, out=vmax)
                    denom = np.sqrt(vmax) + eps
                else:
                    denom = np.sqrt(vh) + eps
                x -= lr * mh / denom
            if t % rec_every == 0:
                rec_t.append(t); rec_x.append(x.copy())
    return np.array(rec_t), np.array(rec_x)

def slope(rec_t, rec_x, schedule="const"):
    """Per-seed OLS slope of x against t (const) or sqrt(t) (decaying), on the
    second half of the run."""
    h = len(rec_t) // 2
    tt = rec_t[h:].astype(float)
    if schedule != "const":
        tt = np.sqrt(tt)
    X = rec_x[h:]
    tc = tt - tt.mean()
    return (tc[:, None] * (X - X.mean(axis=0))).sum(axis=0) / (tc ** 2).sum()

# ----------------------------------------------------------------------------
# RUN
# ----------------------------------------------------------------------------
DELTA = 0.1
DEF = dict(b1=0.9, b2=0.999, eps=1e-8)

log("=== Part 0: theorem bounds at headline point C=1e4 ===")
for C in [1e3, 3e3, 1e4, 1e5]:
    tb = thm_bounds(C, DELTA, DEF["b1"], DEF["b2"])
    log(f"C={C:8.0f}  p={tb['p']:.4e}  J={tb['J']}  (1-p)^J={tb['quiet']:.4f}  "
        f"P(C)={tb['PC']:.5f}  eps_div<={tb['eps_div']:.2f}  eps_conv>={tb['eps_conv']:.1f}")
    OUT.setdefault("thm_bounds", {})[f"C={C:.0f}"] = tb

log("=== Part A: certified drift at defaults (b1=.9, b2=.999, eps=1e-8) ===")
OUT["drift_defaults"] = {}
for C, N in [(1e3, 4_000_000), (3e3, 4_000_000), (1e4, 4_000_000), (1e5, 4_000_000)]:
    r = certified_drift(C, DELTA, DEF["b1"], DEF["b2"], DEF["eps"], N, seed=12345)
    log(f"C={C:8.0f}  D={r['D']:+.5f}  se={r['se']:.5f}  "
        f"bracket=[{r['bracket_lo']:+.5f},{r['bracket_hi']:+.5f}]  "
        f"certified(1-1e-6)=[{r['cert_lo']:+.5f},{r['cert_hi']:+.5f}]  K={r['K']}")
    OUT["drift_defaults"][f"C={C:.0f}"] = r

log("=== Part A2: drift vs eps_tilde sweep at C=1e4 (defaults betas) ===")
C = 1e4
W1, W2, t1, t2, K, p = sample_W(C, DELTA, DEF["b1"], DEF["b2"], 8_000_000, seed=777)
eps_grid = [0.0, 1e-8, 1e-4, 1e-2, 1.0, 3.0, 10.0, 30.0, 100.0, 200.0, 300.0,
            1000.0, 3000.0, 1e4, 3e4, 1e5, 3e5, 1e6]
sweep = []
for eps in eps_grid:
    est, se, lo, hi = drift_from_W(W1, W2, t1, t2, C, DEF["b1"], DEF["b2"], eps)
    sweep.append(dict(eps=eps, D=est, se=se, lo=lo, hi=hi))
    log(f"eps~={eps:10.2e}  D={est:+.6f}  (se {se:.6f})")
OUT["sweep_C1e4"] = sweep
sgn = [s["D"] for s in sweep]
crossings = sum(1 for i in range(len(sgn) - 1) if sgn[i] < 0 <= sgn[i + 1] or sgn[i] >= 0 > sgn[i + 1])
log(f"sign changes on grid: {crossings}")
OUT["sweep_C1e4_crossings"] = crossings
save()

log("=== Part C: phase boundary eps*(C) by bisection, 3 configs ===")
OUT["boundary"] = {}
CONFIGS = [("adam_defaults", 0.9, 0.999), ("adam_llm", 0.9, 0.95), ("rmsprop", 0.0, 0.999)]
for name, b1, b2 in CONFIGS:
    OUT["boundary"][name] = {}
    for C in [1e2, 1e3, 1e4, 1e5]:
        N = 2_000_000 if C <= 1e3 else 6_000_000
        W1, W2, t1, t2, K, p = sample_W(C, DELTA, b1, b2, N, seed=int(1000 + C))
        def D_of(eps):
            return drift_from_W(W1, W2, t1, t2, C, b1, b2, eps)[0]
        d0 = D_of(0.0)
        if d0 >= 0:
            log(f"{name:14s} C={C:8.0f}  D(0)={d0:+.5f} >= 0 : no divergence anywhere")
            OUT["boundary"][name][f"C={C:.0f}"] = dict(eps_star=None, D0=d0)
            continue
        lo_e, hi_e = 1e-6, 1e7
        assert D_of(hi_e) > 0
        for _ in range(50):
            mid = math.sqrt(lo_e * hi_e)
            if D_of(mid) < 0: lo_e = mid
            else: hi_e = mid
        eps_star = math.sqrt(lo_e * hi_e)
        # SE at the crossing -> CI on eps* via local slope
        h = eps_star * 0.2
        dD = (D_of(eps_star + h) - D_of(eps_star - h)) / (2 * h)
        _, se, _, _ = drift_from_W(W1, W2, t1, t2, C, b1, b2, eps_star)
        eps_ci = 1.96 * se / abs(dD) if dD != 0 else float("nan")
        log(f"{name:14s} C={C:8.0f}  D(0)={d0:+.5f}  eps*={eps_star:10.3f}  (+-{eps_ci:.3f})")
        OUT["boundary"][name][f"C={C:.0f}"] = dict(eps_star=eps_star, ci=eps_ci, D0=d0, N=N)
    save()

log("=== Part B: trajectories at C=1e4 (T=2e6, 8 seeds) ===")
C = 1e4; T = 2_000_000; NS = 8; ALPHA = 1e-3
mu = DELTA
runs = [
    ("adam_defaults_const", "adam", DEF["b1"], DEF["b2"], 1e-8, "const", 1.0, T),
    ("adam_defaults_sqrt",  "adam", DEF["b1"], DEF["b2"], 1e-8, "sqrt",  1.0, T),
    ("adam_eps300",         "adam", DEF["b1"], DEF["b2"], 300.0, "const", 1.0, T),
    ("adam_eps3000",        "adam", DEF["b1"], DEF["b2"], 3000.0, "const", 1.0, T),
    ("rmsprop_defaults",    "rmsprop", 0.0, DEF["b2"], 1e-8, "const", 1.0, T),
    ("amsgrad_defaults",    "amsgrad", DEF["b1"], DEF["b2"], 1e-8, "const", 1.0, T),
    ("sgd",                 "sgd", 0, 0, 0, "const", 1.0, T),
    # Loss-scale flip at eps=1e-2 (an RL-style eps), only lambda varies.
    # Boundary eps* ~ 2575 at C=1e4 => flip predicted between lam=1e-5 and 1e-6.
    # Longer T for the two small-|D| runs (statistical power).
    ("adam_eps1e-2_lam1",    "adam", DEF["b1"], DEF["b2"], 1e-2, "const", 1.0, T),
    ("adam_eps1e-2_lam1e-5", "adam", DEF["b1"], DEF["b2"], 1e-2, "const", 1e-5, 2 * T),
    ("adam_eps1e-2_lam1e-6", "adam", DEF["b1"], DEF["b2"], 1e-2, "const", 1e-6, 2 * T),
]
OUT["traj"] = {}
for name, opt, b1, b2, eps, sched, lam, Trun in runs:
    rt, rx = run_traj(opt, C, DELTA, b1, b2, eps, ALPHA, Trun, NS, seed=4242, schedule=sched, lam=lam)
    sl = slope(rt, rx, sched)
    xT = rx[-1]
    log(f"{name:22s} xT: mean={xT.mean():+10.3f} sd={xT.std():8.3f} | "
        f"slope: mean={sl.mean():+.3e} sd={sl.std():.2e}")
    OUT["traj"][name] = dict(xT_mean=float(xT.mean()), xT_sd=float(xT.std()),
                             slope_mean=float(sl.mean()), slope_sd=float(sl.std()),
                             T=Trun, alpha=ALPHA, eps=eps, schedule=sched, lam=lam)
    save()

# theory cross-check: slope should equal -alpha*D (const) resp. -2*alpha*D (sqrt)
d_def = OUT["drift_defaults"]["C=10000"]["D"]
log(f"cross-check: -alpha*D = {-ALPHA*d_def:+.3e}  vs measured "
    f"{OUT['traj']['adam_defaults_const']['slope_mean']:+.3e}")
log(f"cross-check sqrt: -2*alpha*D = {-2*ALPHA*d_def:+.3e}  vs measured "
    f"{OUT['traj']['adam_defaults_sqrt']['slope_mean']:+.3e}")
log(f"cross-check SGD: -alpha*mu = {-ALPHA*mu:+.3e}  vs measured "
    f"{OUT['traj']['sgd']['slope_mean']:+.3e}")

with open("results.json", "w") as f:
    json.dump(OUT, f, indent=1, default=float)
log("DONE -> results.json")