All missions

Motivation

The scalar Bernstein and Bennett inequalities bound the upper tail of a sum of independent, zero-mean random variables that are bounded above: the tail is governed by the total variance at moderate deviations and by the uniform bound in the far tail. Many questions in randomized numerical linear algebra, matrix completion, covariance estimation, graph sparsification and statistical learning reduce to the same question for a sum of independent random matrices, where the quantity of interest is the largest eigenvalue or the spectral norm of the sum rather than a scalar.

Ahlswede and Winter (2002) introduced a matrix Laplace transform method based on the Golden–Thompson inequality. Oliveira (2010) and Gross (2011) proved matrix Bernstein-type bounds with this method. Tropp (arXiv:1004.4389, Found. Comput. Math. 12 (2012)) replaced the Golden–Thompson step with Lieb's concavity theorem, obtaining a master tail bound for independent sums (Theorem 3.6) from which matrix Bennett and Bernstein inequalities follow with the variance parameter $\|\sum_k \mathbb E X_k^2\|$, the natural analogue of the scalar variance. This mission formalizes §6 of that paper.

Timeline:

• 1924–1962, Bernstein and Bennett: scalar inequalities for bounded and subexponential sums.
• 2002, Ahlswede–Winter: the matrix Laplace transform method, via Golden–Thompson.
• 2010, Oliveira; Gross: matrix Bernstein-type bounds by the Ahlswede–Winter method.
• 2010–2012, Tropp: master tail bound via Lieb's theorem; matrix Bennett and Bernstein inequalities (Theorems 6.1, 6.2) with variance parameter $\|\sum_k\mathbb E X_k^2\|$.

Setting

All matrices are $d\times d$ complex matrices with $d\ge 1$. For a self-adjoint (Hermitian) matrix $A$, $\lambda_{\max}(A)$ is its largest eigenvalue and $\|A\|$ its spectral norm. The semidefinite order $A\preccurlyeq B$ means that $B-A$ is positive semidefinite. Functions of a self-adjoint matrix are defined spectrally: if $A=Q\Lambda Q^*$, then $f(A)=Qf(\Lambda)Q^*$; this gives the matrix exponential $e^A$ and, on positive-definite matrices, the matrix logarithm $\log A$. The trace is $\operatorname{tr}$.

A random matrix is a measurable map $X:\Omega\to\mathbb C^{d\times d}$ on a probability space $(\Omega,\mathcal F,\mathbb P)$; its expectation $\mathbb E X$ is taken entrywise. A finite sequence $X_1,\dots,X_n$ of random self-adjoint matrices is independent if the maps are mutually independent.

The setting of the goal: $X_1,\dots,X_n$ are independent random self-adjoint matrices with

$$\mathbb E X_k=0\qquad\text{and}\qquad \lambda_{\max}(X_k)\le R\ \text{almost surely},$$

and the total variance is $\sigma^2:=\big\|\sum_k \mathbb E(X_k^2)\big\|$. Bennett's function is $h(u):=(1+u)\log(1+u)-u$ for $u\ge0$.

Formalization targets

Goal: Theorem 6.1 (matrix Bernstein, bounded case)

For every $t\ge0$, the following chain holds:

$$\mathbb P\Big\{\lambda_{\max}\Big(\sum_k X_k\Big)\ge t\Big\}\le d\cdot\exp\Big(-\frac{\sigma^2}{R^2}\,h\Big(\frac{Rt}{\sigma^2}\Big)\Big)\le d\cdot\exp\Big(\frac{-t^2/2}{\sigma^2+Rt/3}\Big),$$

and the last term is at most $d\,e^{-3t^2/8\sigma^2}$ for $t\le\sigma^2/R$ and at most $d\,e^{-3t/8R}$ for $t\ge\sigma^2/R$. The first inequality is the matrix Bennett inequality, the second the matrix Bernstein inequality, and the case split the split Bernstein inequality.

Milestones, in attack order

1. Theorem 3.6, the master tail bound $\mathbb P\{\lambda_{\max}(\sum_k X_k)\ge t\}\le e^{-\theta t}\operatorname{tr}\exp(\sum_k\log\mathbb E e^{\theta X_k})$ for $\theta>0$.
2. Corollary 3.7: a semidefinite bound $\mathbb E e^{\theta X_k}\preccurlyeq e^{g(\theta)A_k}$ gives $\mathbb P\{\lambda_{\max}(\sum_k X_k)\ge t\}\le d\,e^{-\theta t+g(\theta)\lambda_{\max}(\sum_k A_k)}$.
3. Lemma 6.7: if $\mathbb E X=0$ and $\lambda_{\max}(X)\le1$, then $\mathbb E e^{\theta X}\preccurlyeq\exp\big((e^\theta-\theta-1)\,\mathbb E X^2\big)$ for $\theta>0$.
4. Theorem 6.1(i), the matrix Bennett inequality on its own.
5. The numerical bound $h(u)\ge\dfrac{u^2/2}{1+u/3}$ for $u\ge0$, from the proof of Theorem 6.1.

Further items

Lemma 6.8 (the matrix mgf bound $\mathbb E e^{\theta X}\preccurlyeq\exp(\frac{\theta^2}{2(1-\theta)}A^2)$ under $\mathbb E X^p\preccurlyeq\frac{p!}{2}A^2$) and Theorem 6.2 (matrix Bernstein, subexponential case: $\mathbb P\{\lambda_{\max}(\sum_k X_k)\ge t\}\le d\exp(\frac{-t^2/2}{\sigma^2+Rt})$, split as $d\,e^{-t^2/4\sigma^2}$ and $d\,e^{-t/4R}$) are included as further theorems.

Significance

Matrix Bernstein is the most used inequality of the paper: it bounds the spectral norm of a sum of independent bounded random matrices with only a logarithmic dimensional factor, and it is the standard tool behind sampling bounds for matrix completion, randomized low-rank approximation, spectral sparsification, covariance estimation and the analysis of random features. Theorem 6.1 assumes only an upper bound on the largest eigenvalue of each summand, so it applies to summands with unbounded smallest eigenvalues; applied to $\{X_k\}$ and $\{-X_k\}$ separately it controls both tails.

The results are proved in the paper. As far as the local catalog shows, no machine-checked proof of the matrix Bennett or Bernstein inequality for complex Hermitian matrices exists; the platform has a matrix Bernstein statement for real symmetric matrices in a different form (Wainwright, Theorem 6.17), still open, and scalar Bennett and Bernstein inequalities. Formalizing Theorem 6.1 produces a reusable matrix concentration inequality together with its two ingredients, the matrix mgf bound of Lemma 6.7 and the scalar inequality relating Bennett's and Bernstein's exponents.

Difficulty

The scalar argument multiplies moment generating functions of independent summands. For matrices, $e^{A+B}\ne e^Ae^B$ unless $A$ and $B$ commute, so $\mathbb E\,\operatorname{tr} e^{\theta\sum_k X_k}$ does not factor; the master bound of Theorem 3.6 rests on Lieb's concavity theorem for $A\mapsto\operatorname{tr}\exp(H+\log A)$, which is not in Mathlib. The mgf bound of Lemma 6.7 needs a semidefinite comparison between functions of a random matrix, the transfer rule $f\le g$ on the spectrum $\Rightarrow f(A)\preccurlyeq g(A)$, and the operator monotonicity of the logarithm. The remaining steps (optimizing over $\theta$, the inequality between $h$ and the Bernstein exponent) are real analysis.

Formalization scope

Matrices are Matrix (Fin d) (Fin d) ℂ with [NeZero d]; $\lambda_{\max}$ is the supremum of the real spectrum; $\|\cdot\|$ is the operator norm on $\ell_2^d$; $e^A$ and $\log A$ are Mathlib's continuous functional calculus cfc; the semidefinite order is Mathlib's MatrixOrder; expectation is entrywise. A random matrix is a measurable map, each $X_k$ is Hermitian at every $\omega$, independence is iIndepFun, and the measure is a probability measure. The bound $\lambda_{\max}(X_k)\le R$ is assumed almost surely; no lower bound on the eigenvalues is assumed. The paper's standing regularity (§2.2) is made explicit as integrability of the entries of $X_k$ and $X_k^2$ (and of $e^{\theta X_k}$ where its expectation is taken). The infimum over $\theta>0$ in Theorem 3.6 and Corollary 3.7 is stated for every admissible $\theta>0$.

The printed formulas divide by $R$ and $\sigma^2$, so $0<R$ and $0<\sigma^2$ are hypotheses of Theorems 6.1 and 6.2; with Lean's convention $x/0=0$, the case $\sigma^2=0$ would turn the Bennett bound into the trivial value $d$, so the statement is not allowed to degenerate there. The chain is a conjunction of the probability bound and the inequalities between consecutive bounds, so the weaker reading "the probability is below each bound" does not satisfy it.

A complete development needs the matrix Laplace transform method, Lieb's theorem (or a substitute), operator monotonicity of the logarithm, the transfer rule for the functional calculus, and the scalar inequality $h(u)\ge(u^2/2)/(1+u/3)$. The first four are reusable for every matrix concentration inequality of this paper and beyond. Contributions to any milestone are welcome.

Selected references

• J. A. Tropp, User-Friendly Tail Bounds for Sums of Random Matrices, arXiv:1004.4389v7, 2011; Found. Comput. Math. 12 (2012) 389–434. https://arxiv.org/abs/1004.4389, https://doi.org/10.1007/s10208-011-9099-z
• R. Ahlswede and A. Winter, Strong converse for identification via quantum channels, IEEE Trans. Inform. Theory 48 (2002) 569–579. https://doi.org/10.1109/18.985947
• R. I. Oliveira, Sums of random Hermitian matrices and an inequality by Rudelson, Electron. Commun. Probab. 15 (2010) 203–212. https://arxiv.org/abs/1004.3821
• D. Gross, Recovering low-rank matrices from few coefficients in any basis, IEEE Trans. Inform. Theory 57 (2011) 1548–1566. https://arxiv.org/abs/0910.1879
• E. H. Lieb, Convex trace functions and the Wigner–Yanase–Dyson conjecture, Adv. Math. 11 (1973) 267–288. https://doi.org/10.1016/0001-8708(73)90011-X
• G. Lugosi, Concentration-of-measure inequalities, lecture notes, 2009 (Theorems 5 and 6, scalar Bennett and Bernstein). https://www.econ.upf.edu/~lugosi/anu.pdf

Motivation

Regularized loss minimization over linear predictors is the training problem behind support vector machines, logistic regression, ridge regression and their many variants. Given examples $x_1,\dots,x_n\in\mathbb R^d$, one minimizes an average of convex losses of the predictions $w^\top x_i$ plus a quadratic penalty. Dual coordinate ascent (DCA) solves the Fenchel dual of this problem by optimizing one dual variable at a time; it has been the workhorse of linear SVM solvers such as LIBLINEAR (Hsieh et al., ICML 2008). Before Shalev-Shwartz and Zhang's work, the available analyses of DCA either gave linear rates with constants depending on unknown problem quantities or did not control the duality gap, which is the quantity a practitioner can actually monitor as a stopping criterion.

Shalev-Shwartz and Zhang (2013) analysed stochastic dual coordinate ascent (SDCA), where the coordinate is picked uniformly at random, and proved explicit bounds on the expected duality gap. For Lipschitz losses (hinge loss, absolute loss) the bound is $\tilde O(n+L^2/(\lambda\epsilon))$ iterations, which matches stochastic gradient descent when the accuracy is low and improves on it when $n$ is small compared with $L^2/(\lambda\epsilon)$. This mission formalizes that result, Theorem 1 of the paper, together with the lemmas of its proof.

Timeline. Hsieh, Chang, Lin, Keerthi and Sundararajan (2008) proposed dual coordinate descent for linear SVMs and proved linear convergence with an unspecified constant. Shalev-Shwartz and Zhang (arXiv 2012, JMLR 2013) gave the first explicit duality-gap rates for SDCA, for Lipschitz and for smooth losses. Proximal and accelerated variants followed (Shalev-Shwartz and Zhang, Math. Program. 2016).

Setting

Fix $n\ge1$ examples $x_1,\dots,x_n\in\mathbb R^d$ (Euclidean norm), convex losses $\phi_1,\dots,\phi_n:\mathbb R\to\mathbb R$ and a regularization parameter $\lambda>0$. The primal problem (1) is to minimize

$$P(w)=\frac1n\sum_{i=1}^n\phi_i(w^\top x_i)+\frac\lambda2\|w\|^2 .$$

The convex conjugate of a scalar function is $\phi^*(u)=\sup_z(zu-\phi(z))\in(-\infty,+\infty]$. The dual problem (2) is to maximize, over $\alpha\in\mathbb R^n$,

$$D(\alpha)=\frac1n\sum_{i=1}^n-\phi_i^*(-\alpha_i)-\frac\lambda2\|w(\alpha)\|^2,\qquad w(\alpha)=\frac1{\lambda n}\sum_{i=1}^n\alpha_ix_i .$$

Weak duality gives $D(\alpha)\le P(w)$ for all $\alpha,w$, so the duality gap $P(w(\alpha))-D(\alpha)$ bounds the primal sub-optimality of $w(\alpha)$.

A loss $\phi$ is $L$-Lipschitz if $|\phi(a)-\phi(b)|\le L|a-b|$ for all $a,b$. The paper's standing assumptions are $\|x_i\|\le1$, $\phi_i\ge0$ and $\phi_i(0)\le1$.

Procedure SDCA starts from $\alpha^{(0)}=0$. At iteration $t=1,\dots,T$ it picks a coordinate $i$ uniformly at random, chooses $\Delta\alpha_i$ maximizing

$$-\phi_i^*\bigl(-(\alpha_i^{(t-1)}+\Delta\alpha_i)\bigr)-\frac{\lambda n}{2}\bigl\|w^{(t-1)}+(\lambda n)^{-1}\Delta\alpha_i x_i\bigr\|^2,$$

and sets $\alpha^{(t)}=\alpha^{(t-1)}+\Delta\alpha_ie_i$ and $w^{(t)}=w(\alpha^{(t)})$. The averaging option outputs $\bar\alpha=\frac1{T-T_0}\sum_{t=T_0+1}^T\alpha^{(t-1)}$ and $\bar w=w(\bar\alpha)$. Let $\alpha^*$ maximize $D$ and write $\epsilon_D^{(t)}=D(\alpha^*)-D(\alpha^{(t)})$ for the dual sub-optimality.

Formalization targets

Goal: Theorem 1 (p. 5)

For convex $L$-Lipschitz losses under the standing assumptions, if $\epsilon_P>0$ and

$$T\ \ge\ T_0+n+\frac{4L^2}{\lambda\epsilon_P}\ \ge\ \max\bigl(0,\lceil n\log(0.5\lambda nL^{-2})\rceil\bigr)+n+\frac{20L^2}{\lambda\epsilon_P},$$

then

$$\mathbb E\bigl[P(\bar w)-D(\bar\alpha)\bigr]\le\epsilon_P,\qquad\text{and}\qquad\mathbb E\bigl[D(\alpha^*)-D(\alpha^{(t)})\bigr]\le\epsilon_P/2\ \text{ for all }t\ge T_0 .$$

The constants $4$ and $20$ are the paper's.

Milestones

In the order the proof uses them: the one-coordinate increase (10); the identity $P(w(\alpha))-D(\alpha)=\frac1n\sum_i(\phi_i(w^\top x_i)+\phi_i^*(-\alpha_i)+\alpha_iw^\top x_i)$; Lemma 1 (expected dual increase at least $\frac sn$ times the gap minus $(\frac sn)^2\frac{G}{2\lambda}$); Lemma 2 ($D(\alpha)\le P(w^*)\le P(0)\le1$, $D(0)\ge0$); Lemma 3 ($\phi^*=+\infty$ outside $[-L,L]$); Lemma 4 ($G\le4L^2$); the recursion after (13); the rate (14), $\mathbb E[\epsilon_D^{(t)}]\le\frac{2G}{\lambda(2n+t-t_0)}$; and the averaged duality-gap bound of p. 16.

Significance

Theorem 1 makes the duality gap of SDCA a certified stopping criterion with an explicit iteration count. Because the gap upper-bounds the primal sub-optimality, it also gives a primal guarantee for the averaged output $\bar w$. It covers non-smooth losses, the hinge loss in particular, where linear rates are not available. The analysis pattern, a lower bound on the expected dual increase by the duality gap (Lemma 1) followed by a recursion for the dual sub-optimality, was reused in the analyses of proximal, accelerated and mini-batch dual methods.

The result is proved in the paper; to the best of current knowledge it has no machine-checked proof. The mission produces a formal account of the primal–dual pair of regularized loss minimization with extended-real conjugates, a formal SDCA run with uniformly random coordinates, and the complete chain of lemmas. Lemma 1 is stated with a general strong-convexity modulus $\gamma\ge0$, as in the paper, so it serves the smooth-loss analysis as well.

Difficulty

The obvious first idea is to treat SDCA as stochastic coordinate ascent on a concave function and import a standard rate. That fails: $D$ is neither smooth nor strongly concave for Lipschitz losses, and its domain $\{\alpha:\phi_i^*(-\alpha_i)<\infty\}$ is a box, so the per-coordinate progress cannot be bounded by a gradient norm. The argument instead compares the exact coordinate maximizer with an explicit interpolation towards a sub-gradient, which needs the conjugate's value at a sub-gradient (Fenchel–Young equality) and its infinite values outside $[-L,L]$.

The rate (14) is not a plain geometric recursion. It needs a burn-in phase and then an induction with a step size $s$ that changes with $t$, and the duality gap of the output has to be recovered from the dual progress by Jensen's inequality over the averaged iterates. Everything happens in expectation over the coordinate sequence, while the dual values have to be kept finite along every run.

Formalization scope

Lean conventions:

• $\mathbb R^d$ is EuclideanSpace ℝ (Fin d), examples and losses are indexed by Fin n, and $w^\top x$ is inner ℝ w x.
• The conjugate is an EReal supremum, and $D$ is EReal-valued. Real values of $D$ are used only where the same statement asserts that they are finite.
• The coordinate step is a free map Δ : (Fin n → ℝ) → Fin n → ℝ constrained by an arg-max predicate (IsSDCAStep). Any maximizer is admissible, and one exists for convex real losses.
• The run is a function of the coordinate sequence js : Fin T → Fin n, started at $\alpha^{(0)}=0$, with $w^{(t)}=w(\alpha^{(t)})$.
• Expectations are uniform averages over coordinate sequences (SAGA.Convex.expectIdx, a published definition).
• $\alpha^*$ is any maximizer of $D$, and $w^*$ any minimizer of $P$.
• The constant $G=\max_tG^{(t)}$ of the proof is pinned to its bound $4L^2$. Theorem 1's burn-in $\max(0,\lceil n\log(0.5\lambda nL^{-2})\rceil)$ is exactly the value of $t_0$ obtained with $G=4L^2$ and $\epsilon_D^{(0)}\le1$.
• Lemmas 1 and 4 are stated for a fixed dual-feasible state. The paper's history-averaged form follows from this.

A real-valued conjugate, which returns $0$ where the supremum is infinite, or an unconstrained step map would make the statements false or vacuous. The EReal conjugate and the arg-max predicate rule both out.

Infrastructure needed: Fenchel–Young for scalar convex functions, the existence of sub-gradients of convex real functions and their bound by the Lipschitz constant, Jensen's inequality for the convex $P$ and the concave $D$, and the algebra of uniform averages over $\{1,\dots,n\}^t$. The conjugate and sub-differential definitions and the duality-gap identity can be reused beyond this mission. Proofs of any milestone are welcome, as are lemmas on EReal sums that keep the finiteness bookkeeping short.

Selected references

• S. Shalev-Shwartz and T. Zhang, Stochastic Dual Coordinate Ascent Methods for Regularized Loss Minimization, arXiv:1209.1873v2, 2013; Journal of Machine Learning Research 14:567–599, 2013. https://arxiv.org/abs/1209.1873
• C.-J. Hsieh, K.-W. Chang, C.-J. Lin, S. S. Keerthi and S. Sundararajan, A Dual Coordinate Descent Method for Large-scale Linear SVM, ICML 2008. https://doi.org/10.1145/1390156.1390208
• S. Shalev-Shwartz and T. Zhang, Accelerated Proximal Stochastic Dual Coordinate Ascent for Regularized Loss Minimization, Mathematical Programming 155:105–145, 2016. https://doi.org/10.1007/s10107-014-0839-0

Motivation

The spectral norm of a matrix, its largest singular value, is computable in polynomial time to any accuracy. Its analogue for 3-tensors appears in the analysis of multilinear systems, in the best rank-1 approximation of a data array, in quantum information (the geometric measure of entanglement), and in the injective norm of tensor products. Hillar and Lim, Most Tensor Problems Are NP-Hard (J. ACM 60(6), 2013, Article 45; final version arXiv:0911.1393v5), show that, in contrast with matrices, deciding whether a fixed nonzero rational number is a singular value or the spectral norm of a rational 3-tensor is NP-hard (their Theorems 6.5 and 1.10).

The proof is a reduction from the clique number $\omega(G)$ of a graph. It combines three classical ingredients: the theorem of Motzkin and Straus (1965), which expresses $\omega(G)$ as the value of a quadratic program over the simplex; a theorem of Banach (1938) on the spectral norm of symmetric tensors; and an observation from He, Li and Zhang (2010) on maxima of sums of squared bilinear forms. Deciding $\omega(G)\ge l$ is one of Karp's (1972) NP-complete problems.

Setting

A real 3-tensor is an array $\mathcal A=[\![a_{ijk}]\!]\in\mathbb R^{l\times m\times n}$ with trilinear form $\mathcal A(\mathbf x,\mathbf y,\mathbf z)=\sum_{i,j,k}a_{ijk}x_iy_jz_k$. With the Euclidean norm $\|\cdot\|_2$, its spectral norm is

$$\|\mathcal A\|_{2,2,2}=\sup_{\mathbf x,\mathbf y,\mathbf z\ne\mathbf 0}\frac{|\mathcal A(\mathbf x,\mathbf y,\mathbf z)|}{\|\mathbf x\|_2\|\mathbf y\|_2\|\mathbf z\|_2}.$$

A real $\sigma$ is a unit $\ell^2$-singular value of $\mathcal A$ if there are unit vectors $\mathbf u,\mathbf v,\mathbf w$ with

$$\textstyle\sum_{j,k}a_{ijk}v_jw_k=\sigma u_i,\quad \sum_{i,k}a_{ijk}u_iw_k=\sigma v_j,\quad \sum_{i,j}a_{ijk}u_iv_j=\sigma w_k$$

for all $i,j,k$: the stationarity conditions of $\mathcal A(\mathbf u,\mathbf v,\mathbf w)$ on a product of unit spheres.

Let $G=(V,E)$ be a simple graph on $V=\{1,\dots,v\}$, $v\ge1$, with $e$ edges $\{i_k,j_k\}$, $k=1,\dots,e$, and clique number $\omega$. Let $\Delta_v$ be the standard simplex, $A_G$ the adjacency matrix and $J$ the all-ones matrix. For a positive integer $l$ set $Q_l=A_G+\frac1lJ$ and $M_l=\max_{\mathbf x\in\Delta_v}\mathbf x^\top Q_l\mathbf x$. Let $E_k=\frac12E_{i_kj_k}+\frac12E_{j_ki_k}$, and let $T_l$ be the maximum over unit $\mathbf u,\mathbf v\in\mathbb R^v$, $\mathbf w\in\mathbb R^{l+2e}$ of the multilinear form

$$\sum_{i=1}^l\Big(\mathbf u^\top\tfrac1lI\mathbf v\Big)w_i+\sum_{k=1}^e(\mathbf u^\top E_k\mathbf v)w_{l+k}+\sum_{k=1}^e(\mathbf u^\top E_k\mathbf v)w_{e+l+k}. \tag{26}$$

The graph tensor $\mathcal A_l\in\mathbb Q^{v\times v\times(l+2e)}$ has as entry $a_{ijk}$ the coefficient of $u_iv_jw_k$ in (26).

Formalization targets

Goal: Theorems 6.5 and 1.10, encoding-free core

For every graph $G$ on $v\ge1$ vertices and every enumeration of its edges:

$$\forall l\ge1:\ \|\mathcal A_l\|_{2,2,2}=1\iff l=\omega;\qquad 1\text{ is a unit singular value of }\mathcal A_\omega;\qquad \forall l>\omega:\ 1\text{ is not one of }\mathcal A_l,$$

so $\omega$ is the largest $l\in\{1,\dots,v\}$ for which $1$ is a unit singular value of $\mathcal A_l$. The paper's $v$-query Turing reductions rest on exactly this.

Milestones (in attack order)

1. Theorem 6.7 (Motzkin–Straus): $1-1/\omega=2\max_{\mathbf x\in\Delta_v}\sum_{\{i,j\}\in E}x_ix_j$.
2. (22): $M_l=1+\frac{\omega-l}{l\omega}$; so $M_\omega=1$, $M_l>1$ for $l<\omega$, $M_l<1$ for $l>\omega$.
3. Lemma 6.8: $M_l=N_l$, where $N_l=\max_{\|\mathbf u\|_2=1}\{\sum_{i=1}^l(\mathbf u^\top\frac1lI\mathbf u)^2+2\sum_k(\mathbf u^\top E_k\mathbf u)^2\}$.
4. Theorem 6.9 (Banach), (24): for a symmetric 4-tensor, $\|\mathcal S\|_{2,2,2,2}=\sup_{\mathbf x\neq\mathbf 0}|\mathcal S(\mathbf x,\mathbf x,\mathbf x,\mathbf x)|/\|\mathbf x\|_2^4$.
5. Proposition 6.10 (He–Li–Zhang): for symmetric $A_k$, $\max_{\|\mathbf u\|=\|\mathbf v\|=1}\sum_k(\mathbf u^\top A_k\mathbf v)^2=\max_{\|\mathbf v\|=1}\sum_k(\mathbf v^\top A_k\mathbf v)^2$.
6. Lemma 6.11: $T_l=M_l^{1/2}$, hence $T_l=1$ iff $l=\omega$, $T_l>1$ iff $l<\omega$, $T_l<1$ iff $l>\omega$.
7. Proof of Theorem 6.5: $\|\mathcal A_l\|_{2,2,2}=T_l$, $T_l$ is a unit singular value of $\mathcal A_l$, and every unit singular value has absolute value at most $T_l$.
8. §7 (content of Theorem 1.13): $\min_{\sigma\ge0,\ \text{unit }\mathbf u,\mathbf v,\mathbf w}\|\mathcal A-\sigma\,\mathbf u\otimes\mathbf v\otimes\mathbf w\|_F^2=\|\mathcal A\|_F^2-\|\mathcal A\|_{2,2,2}^2$, attained only with $\sigma=\|\mathcal A\|_{2,2,2}$.

Significance

The result. The goal shows that a single comparison with $1$ of the spectral norm, or of the singular-value spectrum, of an explicit rational tensor of size $v\times v\times(l+2e)$ encodes whether $l=\omega(G)$. It is the step that makes the spectral norm, the largest singular value, and (through §7) the best rank-1 approximation of a 3-tensor NP-hard, whereas all three are polynomial-time for matrices. It explains why the higher-order power method and its relatives carry no global guarantees, and it is the starting point for inapproximability results such as the paper's Theorem 1.11.

Formalizing it. The results are proved on paper; none of them is formalized. Motzkin–Straus, Banach's polarization theorem for symmetric 4-tensors, and the He–Li–Zhang identity are classical results that are not in Mathlib and not on the platform; each is reusable well beyond this mission (Turán-type bounds, polynomial optimization over spheres, injective norms). The mission also checks the paper's construction at the level of indices and constants.

Difficulty

The construction of $\mathcal A_l$ is explicit, and the comparison with $1$ is arithmetic once $T_l=M_l^{1/2}$ is known. Everything rests on the three cited theorems. Motzkin–Straus is an optimization statement over the simplex whose maximizers may lie on the boundary and need not be unique, so a first-order computation at an interior critical point does not settle it. Banach's theorem asserts equality, not equality up to a constant: the standard polarization identity bounds the 4-linear form by its diagonal values only up to a factor. Proposition 6.10 is where the paper's own argument is loose: the 4-linear form it introduces is not symmetric, so (24) does not apply as printed. Finally, the singular-value clauses require that the maximum of a trilinear form on a product of spheres is a stationary value. This must be proved, not assumed.

Formalization scope

Tensors are arrays Fin l → Fin m → Fin n → ℝ with 0-based indices: the paper's $a_{ijk}$ is A (i-1) (j-1) (k-1). The Euclidean norm is Real.sqrt (∑ i, x i ^ 2). Suprema and maxima are sSup of sets that are nonempty and bounded under the stated hypotheses; where the paper's "max" carries content (Motzkin–Straus, (22), §7), attainment is stated with IsGreatest/IsLeast. The graph is a SimpleGraph (Fin v) with $v\ge1$, clique number Mathlib's SimpleGraph.cliqueNum, and its edges are enumerated by an arbitrary Fin e ≃ G.edgeSet. Every statement holds for each enumeration, so the tensor depends on the choice and the conclusions do not. $\mathcal A_l$ is a rational tensor, read over $\mathbb R$ by casting.

Explicit readings of the paper:

• Definition 6.1 has no norm constraint. Since $(\sigma,\mathbf u,\mathbf v,\mathbf w)\mapsto(t\sigma,t\mathbf u,t\mathbf v,t\mathbf w)$ preserves (20), unnormalized singular values make every nonzero real a singular value of $\mathcal A_l$, so the clause "not for $l>\omega$" would be false. Singular vectors are therefore unit vectors, as in the paper's derivation on p. 0:21.
• The statement of Lemma 6.11 prints $w_{m+l+k}$. There is no $m$ in §6, and the proof uses $w_{e+l+k}$, which is formalized.
• $E_k$ is written entrywise as $\frac12[\{i,j\}=\{i_k,j_k\}]$, which equals $\frac12E_{i_kj_k}+\frac12E_{j_ki_k}$ for a loopless edge.
• "$T_l$ is just the maximum $\ell^2$-singular value of $\mathcal A_l$" is made explicit as milestone 7.
• The goal does not claim that $1$ fails to be a singular value of $\mathcal A_l$ for $l<\omega$; the paper neither claims nor uses this.

This is the mathematical core of a polynomial-time Turing reduction, not an NP-hardness theorem. The polynomial size of $\mathcal A_l$, the $v$-query procedure, the scaling to a general fixed $0\ne\sigma\in\mathbb Q$, and the NP-completeness of CLIQUE (Karp 1972, cited) are not formalized. Theorem 1.11 and Corollary 1.12 are out of scope.

Two formalizations would make the goal trivial, and both are excluded. One is dropping the unit-norm constraint, which makes every nonzero number a singular value. The other is a spectral norm taken as a supremum over all vectors without normalization, which makes it unbounded and junk-valued; here the quotients are normalized.

Contributions are welcome on each milestone. Motzkin–Straus, Banach (24) and He–Li–Zhang are self-contained and reusable. Lemma 6.11 and milestone 7 need Cauchy–Schwarz and a Lagrange-type stationarity argument on spheres, which Mathlib supports through IsLocalExtrOn and the Lagrange multiplier files.

Selected references

• C. J. Hillar, L.-H. Lim, Most Tensor Problems Are NP-Hard, J. ACM 60(6), Article 45, 2013. https://doi.org/10.1145/2512329 (final version arXiv:0911.1393v5, https://arxiv.org/abs/0911.1393v5)
• T. S. Motzkin, E. G. Straus, Maxima for graphs and a new proof of a theorem of Turán, Canadian J. Math. 17, 533–540, 1965. https://doi.org/10.4153/CJM-1965-053-6
• S. Banach, Über homogene Polynome in (L²), Studia Math. 7, 36–44, 1938. https://doi.org/10.4064/sm-7-1-36-44
• S. He, Z. Li, S. Zhang, Approximation algorithms for homogeneous polynomial optimization with quadratic constraints, Math. Program. Ser. B 125(2), 353–383, 2010. https://doi.org/10.1007/s10107-010-0409-z
• R. M. Karp, Reducibility among combinatorial problems, in Complexity of Computer Computations, Plenum, 85–103, 1972. https://doi.org/10.1007/978-1-4684-2001-2_9

Motivation

Conic optimization over the exponential cone $K_{\exp}=\operatorname{cl}\{x\in\mathbb R^3: x_1\ge x_2e^{x_3/x_2},\ x_2>0\}$ models entropy, relative entropy, logistic regression, log-sum-exp and geometric programming. Unlike the nonnegative orthant, the second-order cone and the semidefinite cone, the exponential cone is not symmetric (self-scaled), so the Nesterov–Todd primal-dual machinery that makes solvers such as SeDuMi and MOSEK fast does not apply as is. Dahl and Andersen (Math. Program. 194 (2022) 341–370) describe the nonsymmetric-cone algorithm implemented in MOSEK (version 9.2 in the paper's experiments). Its central new ingredient is a corrector for nonsymmetric cones, an analogue of Mehrotra's predictor-corrector that is responsible for much of the practical performance of symmetric-cone solvers. The paper reports that the corrector gives a substantial and consistent reduction in iteration counts across its test sets.

Timeline. Nesterov and Todd (1997, 1998) built primal-dual methods for self-scaled cones. Nesterov, Todd and Ye (1999) and Tunçel (2001) studied primal-dual scalings for general cones, Tunçel proving polynomial complexity of an infeasible method without a corrector under bounded scalings. Skajaa and Ye (2015) gave a homogeneous method for nonsymmetric cones with a Runge–Kutta corrector. Myklebust and Tunçel (2014) analysed BFGS-type scalings. Dahl and Andersen (2022) introduced the third-order corrector studied in this mission.

Setting

Let $A\colon\mathbb R^n\to\mathbb R^m$ be linear with full row rank, $b\in\mathbb R^m$, $c\in\mathbb R^n$, and let $\hat K\subseteq\mathbb R^n$ be a proper cone (pointed, closed, convex, with nonempty interior). A $\hat\vartheta$-logarithmically homogeneous self-concordant barrier (LHSCB) for $\hat K$ is a $C^3$ convex function $\hat F$ on $\operatorname{int}\hat K$ with $|\hat F'''(x)[u,u,u]|\le2(\hat F''(x)[u,u])^{3/2}$ and $\hat F(\tau x)=\hat F(x)-\hat\vartheta\log\tau$. The primal problem minimizes $\langle c,\hat x\rangle$ subject to $A\hat x=b$, $\hat x\in\hat K$; the dual maximizes $\langle b,y\rangle$ subject to $c-A^Ty=\hat s\in\hat K^*$.

The homogeneous model adds $\tau,\kappa\ge0$: $x=(\hat x,\tau)$, $s=(\hat s,\kappa)$, $K=\hat K\times\mathbb R_+$, $F(x)=\hat F(\hat x)-\log\tau$, $\vartheta=\hat\vartheta+1$. With $z=(x,s,y)$, the residual is

$$G(z)=\bigl(A\hat x-b\tau,\ -A^Ty+c\tau-\hat s,\ b^Ty-c^T\hat x-\kappa\bigr),$$

a linear map whose $(y,\hat x,\tau)$-block is skew-symmetric. An iterate has $x\in\operatorname{int}K$, $s\in\operatorname{int}K^*$, and $\mu=\langle x,s\rangle/\vartheta$. The shadow iterates are $\tilde x=-F_*'(s)$ and $\tilde s=-F'(x)$, with $F_*$ the conjugate barrier. A scaling is a nonsingular $W$ with $v=Wx=W^{-T}s$ and $\tilde v=W\tilde x=W^{-T}\tilde s$ (the double secant equations (8)).

The affine direction $\Delta z^a$ solves $G(\Delta z^a)=-G(z)$, $W\Delta x^a+W^{-T}\Delta s^a=-v$ (11). The corrector is $\eta=-\tfrac12F'''(x)[\Delta x^a,(F''(x))^{-1}\Delta s^a]$ (16). For a centering parameter $\gamma>0$, the combined direction $\Delta z$ solves

$$G(\Delta z)=-(1-\gamma)G(z),\qquad W\Delta x+W^{-T}\Delta s=-v+\gamma\mu\tilde v-W^{-T}\eta. \tag{18}$$

Formalization targets

Goal: Lemma 4 (p. 353)

Every solution of (18) satisfies

$$\langle s,\Delta x\rangle+\langle x,\Delta s\rangle=-(1-\gamma)\langle x,s\rangle,\qquad\langle\Delta x,\Delta s\rangle=0,$$

and for all $\alpha\in\mathbb R$

$$G(z+\alpha\Delta z)=(1-\alpha(1-\gamma))G(z),\qquad\langle x+\alpha\Delta x,s+\alpha\Delta s\rangle=(1-\alpha(1-\gamma))\langle x,s\rangle.$$

The goal fixes no constants beyond those of the direction itself. It holds for every solution of (18), not a chosen one.

Milestones

1. §2 homogeneity identities (p. 345): $F'(\tau x)=F'(x)/\tau$, $F''(\tau x)=F''(x)/\tau^2$, $F''(x)x=-F'(x)$, $F'''(x)[x]=-2F''(x)$, $\langle F'(x),x\rangle=-\vartheta$.
2. §2 sum rule for the homogenizing block (pp. 345, 347): $\hat F(\hat x)-\log\tau$ is a $(\hat\vartheta+1)$-LHSCB for $\hat K\times\mathbb R_+$.
3. Lemma 2 (p. 351): the affine direction satisfies $\langle s,\Delta x^a\rangle+\langle x,\Delta s^a\rangle=-\langle x,s\rangle$, $\langle\Delta x^a,\Delta s^a\rangle=0$, and $\langle x+\alpha\Delta x^a,s+\alpha\Delta s^a\rangle=(1-\alpha)\langle x,s\rangle$.
4. Lemma 3 (p. 352): the pure corrector direction ($G(\Delta z^c)=0$, $W\Delta x^c+W^{-T}\Delta s^c=-W^{-T}\eta$) satisfies $\langle s,\Delta x^c\rangle+\langle x,\Delta s^c\rangle=0$ and $\langle\Delta x^c,\Delta s^c\rangle=0$.

Significance

The result. Lemma 4 says that a step of any length along the combined direction reduces the infeasibility residuals and the complementarity gap by exactly the same factor. Starting from a point with $\mu^0=1$, the iterates satisfy $G(z^k)=\mu^kG(z^0)$ and $\langle x^k,s^k\rangle=\mu^k\vartheta$, so the algorithm needs no merit function to balance feasibility against optimality, in contrast with earlier nonsymmetric methods. The lemma also shows that adding the third-order corrector does not cost this property, which is what justifies using the corrector in MOSEK's exponential-cone solver.

Formalizing it. The result is proved in the paper; nothing here is open. To our knowledge none of Lemmas 2–4, the §2 homogeneity identities or the barrier sum rule has a machine-checked proof. The mission produces a formal homogeneous model for a general proper cone with a general LHSCB and a formal account of third-derivative calculus for logarithmically homogeneous barriers. Both are reusable in any analysis of interior-point methods for nonsymmetric cones (power cones, generalized power cones, hyperbolicity cones).

Difficulty

The obvious reading of Lemma 4 as bookkeeping with the linear map $G$ stalls at the corrector: the term $W^{-T}\eta$ in (18) involves the third derivative of the barrier, and nothing in linear algebra says how it interacts with the complementarity gap. Controlling it needs the third-derivative calculus of logarithmically homogeneous barriers. In Lean that means relating gradient, the Hessian as fderiv of the gradient, and the third derivative as fderiv of the Hessian, and working with the symmetry of iterated Fréchet derivatives of a $C^3$ function, for which the Mathlib interface is thin. A second obstacle is the augmented barrier: the barrier properties of $\hat F$ must be transported to $\hat F(\hat x)-\log\tau$ on $\mathbb R^{n+1}$ through the splitting $\mathbb R^{n+1}=\mathbb R^n\times\mathbb R$, including the boundary behaviour at the corner $\hat x\in\partial\hat K$, $\tau=0$.

The first identity of Lemma 4 is misprinted in the paper as $=0$. A formalization following the printed statement would be false, not merely awkward.

Formalization scope

Vectors of the augmented space are EuclideanSpace ℝ (Fin (n + 1)); $\hat x$ is the first $n$ coordinates and $\tau$ (resp. $\kappa$) the last. $A$ is a continuous linear map and $A^T$ its Hilbert adjoint; full row rank is surjectivity. Proper cones, the Hessian hess, the barrier class IsLogHomBarrier and the conjugate conj are the published SelfScaledIPM.ShortStep.Setting. That barrier class adds three standard conditions to the paper's two (boundary blow-up, the $\vartheta$-inequality and a positive definite Hessian), which the paper uses implicitly when it writes $(F''(x))^{-1}$ and calls $F$ a barrier. $W$ is a continuous linear equivalence and $W^{-T}$ the adjoint of $W^{-1}$. $F'''(x)[u,w]$ is fderiv ℝ (hess F) x u w and $(F''(x))^{-1}$ is ContinuousLinearMap.inverse. Every inner product includes the homogenizing term $\tau\kappa$.

Repairs, generalizations and specializations, each disclosed in the item statements:

• Repair: Lemma 4's first identity is stated as $-(1-\gamma)\langle x,s\rangle$, as derived in the paper's own proof, instead of the printed $0$.
• Generalization: one proper cone $\hat K$ with one barrier replaces the paper's product $K_1\times\dots\times K_k$. The product with the sum barrier is a special case.
• Generalization: $W$ is any nonsingular map satisfying (8). The paper's block-diagonal scaling with $W_{k+1}=\sqrt{\kappa/\tau}$ is a special case.
• Specialization: the sum rule is stated for $K_2=\mathbb R_+$, $F_2=-\log$, $\vartheta_2=1$, the only case the paper applies.

The directions are predicates ("$\Delta z$ solves (11)"), and the lemmas quantify over all solutions. The hypotheses are satisfiable. For the linear program $n=m=1$, $A=1$, $b=c=1$, $\hat K=\mathbb R_+$, $\hat F=-\log$, $x=s=(1,1)$, $y=0$, $W=I$, system (11) has the unique solution $\Delta x^a=(0,0)$, $\Delta s^a=(-1,-1)$, $\Delta y^a=1$. Statements that hold only because a junk value appears are ruled out: every statement requires $x\in\operatorname{int}K$ and $s\in\operatorname{int}K^*$, where $\log\tau$, the inverse Hessian and the conjugate barrier take their true values. A formalization that drops $\tau,\kappa$ from the inner products, or that fixes one solution of (18) instead of quantifying over all, would be a different lemma and does not count.

Welcome contributions: the general sum rule for products of cones, a Fin.append splitting library for EuclideanSpace, and proofs of the homogeneity identities that can be reused by the series' other missions.

Selected references

• J. Dahl, E. D. Andersen, A primal-dual interior-point algorithm for nonsymmetric exponential-cone optimization, Math. Program. 194 (2022) 341–370. https://doi.org/10.1007/s10107-021-01631-4
• Yu. Nesterov, M. J. Todd, Primal-dual interior-point methods for self-scaled cones, SIAM J. Optim. 8 (1998) 324–364. https://doi.org/10.1137/S1052623495290209
• L. Tunçel, Generalization of primal-dual interior-point methods to convex optimization problems in conic form, Found. Comput. Math. 1 (2001) 229–254. https://doi.org/10.1007/s002080010009
• A. Skajaa, Y. Ye, A homogeneous interior-point algorithm for nonsymmetric convex conic optimization, Math. Program. 150 (2015) 391–422. https://doi.org/10.1007/s10107-014-0773-1
• T. Myklebust, L. Tunçel, Interior-point algorithms for convex optimization based on primal-dual metrics, arXiv:1411.2129 (2014). https://arxiv.org/abs/1411.2129

Why adaptive timesteps

The explicit Euler–Maruyama scheme is the standard way to simulate a stochastic differential equation

$$dX_t=f(X_t)\,dt+g(X_t)\,dW_t,\qquad X_0=x_0,$$

on $\mathbb R^m$ driven by a $d$-dimensional Brownian motion. When $f$ and $g$ are globally Lipschitz it converges strongly with order $\tfrac12$ in the step size (Kloeden and Platen 1992). Many models of interest in physics, chemistry and statistics have a drift that grows faster than linearly but points inwards, such as $f(x)=-x^3$. For these the SDE is well posed and its solution has finite moments of all orders, but the uniform-step scheme does not: for $dX_t=-X_t^3dt+dW_t$, Hutzenthaler, Jentzen and Kloeden proved that $\mathbb E\|X_T\|^p\to\infty$ as the step size $h\to0$, for every $T>0$ and $p\ge2$ (Proc. R. Soc. A 2011).

Several remedies have been proposed: implicit schemes (Higham, Mao and Stuart 2002; Mao and Szpruch 2013), the tamed Euler scheme (Hutzenthaler, Jentzen and Kloeden 2012), truncated schemes (Mao 2015). Fang and Giles (Ann. Appl. Probab. 2020) keep the explicit Euler–Maruyama update and change only the step: the timestep is a function of the current state, small where the drift is large. This mission formalizes their finite-time results: stability of the adaptive scheme (Theorem 1) and strong convergence of order $\tfrac12$ (Theorem 3).

The adaptive scheme

Let $\|x\|$ and $\langle x,y\rangle$ be the Euclidean norm and inner product on $\mathbb R^m$, and $\|A\|=(\sum_{i,j}A_{ij}^2)^{1/2}$ the Frobenius norm of an $m\times d$ matrix. Given a timestep function $h:\mathbb R^m\to(0,\infty)$, the scheme (5) starts from $t_0=0$, $\widehat X_0=x_0$ and sets $h_n=h(\widehat X_{t_n})$,

$$t_{n+1}=t_n+h_n,\qquad \widehat X_{t_{n+1}}=\widehat X_{t_n}+f(\widehat X_{t_n})h_n+g(\widehat X_{t_n})(W_{t_{n+1}}-W_{t_n}).$$

The grid is random: it depends on the path. For a time $t$, $\underline t$ is the last grid time not after $t$; the piecewise constant interpolant is $\overline X_t=\widehat X_{\underline t}$ and the continuous interpolant is

$$\widehat X_t=\widehat X_{\underline t}+f(\widehat X_{\underline t})(t-\underline t)+g(\widehat X_{\underline t})(W_t-W_{\underline t}).$$

The hypotheses (Assumptions 1–4, pp. 528–530):

1. Assumption 1: $f,g$ locally Lipschitz, $\langle x,f(x)\rangle\le\alpha_1\|x\|^2+\beta_1$ and $\|g(x)\|^2\le\alpha_1\|x\|^2+\beta_1$.
2. Assumption 2: $h$ continuous and positive with $\langle x,f(x)\rangle+\tfrac12h(x)\|f(x)\|^2\le\alpha\|x\|^2+\beta$ for some $\alpha,\beta>0$.
3. Assumption 3: a family $h^\delta$, $0<\delta\le1$, with $\delta\min(T,h(x))\le h^\delta(x)\le\min(\delta T,h(x))$; $\delta\to0$ plays the role of $h\to0$.
4. Assumption 4: $\langle x-y,f(x)-f(y)\rangle\le\tfrac12\alpha\|x-y\|^2$, $\|g(x)-g(y)\|^2\le\tfrac12\alpha\|x-y\|^2$, and $\|f(x)-f(y)\|\le(\gamma(\|x\|^q+\|y\|^q)+\mu)\|x-y\|$.

For $f(x)=-x^3$ in one dimension, $h(x)=\min(1,1/x^2)$ satisfies Assumption 2.

Formalization targets

Goal: Theorem 3 (strong convergence order, p. 531)

Under Assumption 4, with $h^\delta$ satisfying Assumption 3 for an $h$ satisfying Assumption 2, for every $p>0$ there is $C_{p,T}$ such that for all $\delta\in(0,1]$

$$\mathbb E\Big[\sup_{0\le t\le T}\|\widehat X_t-X_t\|^p\Big]\le C_{p,T}\,\delta^{p/2}.$$

Milestones

• Lemma 1 (p. 529): under Assumption 1, $\mathbb E[\sup_{0\le t\le T}\|X_t\|^p]<\infty$ for all $p>0$.
• (35), (37) (p. 548): one-step and partial-step energy inequalities for the projected K-scheme $\widehat X^K_{t_{n+1}}=P_K(\widehat X^K_{t_n}+f(\widehat X^K_{t_n})h_n+g(\widehat X^K_{t_n})\Delta W_n)$, $P_K(Y)=\min(1,K/\|Y\|)Y$.
• §6.1 Step 3 (p. 550): $\mathbb E[\sup_{0\le t\le T}\|\widehat X^K_t\|^p]\le C_{p,T}$ uniformly in $K>\|x_0\|$.
• Theorem 1 (p. 529): $T$ is almost surely attained and $\mathbb E[\sup_{0\le t\le T}\|\widehat X_t\|^p]<C_{p,T}$, with $C_{p,T}$ independent of $h$.
• §6.2, after (41) (p. 553): $\mathbb E\|\widehat X_s-\overline X_s\|^{2p}\le C^3_{p,T}\delta^p$.
• (40) (p. 552): the integral inequality for $e_t=\widehat X_t-X_t$ to which Grönwall's inequality is applied.

Significance

Theorem 1 shows that an explicit scheme, at the cost of one drift and one diffusion evaluation per step, is stable for one-sided Lipschitz drifts of polynomial growth, the class where the uniform-step scheme fails. Theorem 3 shows that refining the timestep function keeps the classical order $\tfrac12$. Combined with the paper's bound on the expected number of steps (Lemma 2) this gives order $\tfrac12$ in terms of computational cost, and the results underpin the paper's adaptive multilevel Monte Carlo estimators. The uniform-in-$h$ form of Theorem 1 is what lets Theorem 3 apply the stability bound to the whole family $h^\delta$.

The paper's proofs are complete on paper; none of these results has a machine-checked proof. Formalizing them needs a theory of Euler schemes on a random, state-dependent grid, which is not in Mathlib and is not covered by uniform-grid formalizations of tamed or implicit schemes. The definitions here (the scheme on random times, the projected K-scheme) are reusable for the infinite-time companion mission and for other adaptive schemes.

Difficulty

The obvious argument for the uniform-step scheme bounds moments step by step, using deterministic step sizes and independent Gaussian increments. Here the step $h_n$ depends on $\widehat X_{t_n}$, so the grid times are stopping times and the increments $W_{t_{n+1}}-W_{t_n}$ have random length. Sums over steps become Itô integrals only after the scheme is written as the solution of $d\widehat X_t=f(\overline X_t)dt+g(\overline X_t)dW_t$, and moment bounds for those integrals need the Burkholder–Davis–Gundy inequality on a random grid. A second difficulty is that $T$ need not be reached: if $h(\widehat X_{t_n})$ decays fast enough, $t_n$ could converge below $T$. Excluding this needs a positive lower bound on $h$ over bounded sets together with moment bounds, which the paper obtains by first analysing a version of the scheme projected onto a ball.

Formalization scope

Lean conventions: $\mathbb R^m$ is the published SDEState m (EuclideanSpace ℝ (Fin m)); matrices are the published Diffusion m d, whose norm is the Frobenius norm; $W$ is a standard $d$-dimensional $(\mathcal F_t)$-Brownian motion in the sense of the published IsWienerMartingale; the SDE solution is any process satisfying the published IsSolution on $[0,T]$ (existence and uniqueness are not claimed). Times are in $[0,\infty)$, indices start at $0$, and $x_0$ is deterministic. The scheme is defined pathwise; the index of $\underline t$ is a junk value $0$ when the grid never passes $t$, an event of probability zero (Theorem 1). $h$ and $h^\delta$ are real-valued and their positivity is a hypothesis; each $h^\delta$ is assumed measurable. Moments are $[0,\infty]$-valued integrals of extended norms, never Bochner integrals. Every existential constant is chosen before the parameter it must not depend on: before $\delta$ in Theorem 3 and its milestones, before $h$ (and $K$) in Theorem 1 and Step 3. Results the paper proves for $p\ge4$ inside its proofs are stated for $p\ge4$. Pages are the journal's printed pages (PDF page + 525).

A formalization that places the constant after $\delta$, uses a deterministic grid, or allows $h=0$ would be trivially true or a different theorem; these are ruled out by the statements.

Contributions welcome: proofs of the milestones, and of general facts they need (Burkholder–Davis–Gundy, moments of Brownian increments over stopping-time intervals, Grönwall's inequality for $[0,\infty]$-valued functions).

Selected references

• W. Fang, M. B. Giles, Adaptive Euler–Maruyama method for SDEs with nonglobally Lipschitz drift, Ann. Appl. Probab. 30(2), 526–560, 2020. https://doi.org/10.1214/19-AAP1507 (preprint arXiv:1609.08101)
• D. J. Higham, X. Mao, A. M. Stuart, Strong convergence of Euler-type methods for nonlinear stochastic differential equations, SIAM J. Numer. Anal. 40, 1041–1063, 2002. https://doi.org/10.1137/S0036142901389530
• M. Hutzenthaler, A. Jentzen, P. E. Kloeden, Strong and weak divergence in finite time of Euler's method for stochastic differential equations with non-globally Lipschitz continuous coefficients, Proc. R. Soc. A 467, 1563–1576, 2011. https://doi.org/10.1098/rspa.2010.0348
• M. Hutzenthaler, A. Jentzen, P. E. Kloeden, Strong convergence of an explicit numerical method for SDEs with nonglobally Lipschitz continuous coefficients, Ann. Appl. Probab. 22, 1611–1641, 2012. https://doi.org/10.1214/11-AAP803
• X. Mao, L. Szpruch, Strong convergence and stability of implicit numerical methods for stochastic differential equations with non-globally Lipschitz continuous coefficients, J. Comput. Appl. Math. 238, 14–28, 2013. https://doi.org/10.1016/j.cam.2012.08.015
• X. Mao, The truncated Euler–Maruyama method for stochastic differential equations, J. Comput. Appl. Math. 290, 370–384, 2015. https://doi.org/10.1016/j.cam.2015.06.002
• P. E. Kloeden, E. Platen, Numerical Solution of Stochastic Differential Equations, Springer, 1992. https://doi.org/10.1007/978-3-662-12616-5