Numerical Analysis

Motivation

Many-server queues are the basic model of call centres and hospital wards, and the Halfin–Whitt regime (the number of servers $n$ grows while the traffic intensity approaches one at rate $n^{-1/2}$) is where such systems operate efficiently and with short delays. In this regime the centred and scaled queue process of a $G/Ph/n+GI$ queue (general arrivals, phase-type service, $n$ servers, generally distributed patience) converges to a $d$-dimensional diffusion, $d$ being the number of service phases (Dai, He and Tezcan, Ann. Appl. Probab. 2010). For $GI/Ph/n+M$ queues the stationary distribution of the queue converges to the invariant measure $\mu$ of that diffusion (Dai, Dieker and Gao, Queueing Systems 2014). So $\mu$ is the object one would like to compute.

With exponential service the diffusion is one-dimensional and $\mu$ is explicit. With phase-type service it is multidimensional with a piecewise-linear drift; Dieker and Gao (Ann. Appl. Probab. 2013) proved that it is exponentially ergodic with a unique $\mu$, but no closed form is known, and the characterization of multidimensional piecewise-linear diffusions was left open in Browne and Whitt (1995). Jin, Pang, Xu and Xu (arXiv:2209.07361) approximate $\mu$ by the invariant measure of the Euler–Maruyama discretization and give an explicit error rate in the Wasserstein-1 distance.

Setting

Fix $d\ge1$, a patience rate $\alpha>0$, a slack $\beta\in\mathbb R$, a probability vector $p\in\mathbb R^d$, service rates $v_1,\dots,v_d>0$ and a sub-stochastic phase-transition matrix $P$ with $P_{ii}=0$ and $I-P$ invertible. Let $R=(I-P')\,\mathrm{diag}(v)$, $e=(1,\dots,1)'$, and assume the phase-type distribution has mean one, $e'R^{-1}p=1$. The drift is the piecewise-linear map

$$g(x)=-\beta p-Rx+(R-\alpha I)p\,(e'x)^+,\qquad y^+=\max\{0,y\}.$$

Let $\sigma$ be a $d\times d$ matrix with $\xi'\sigma\sigma'\xi\ge c\,\xi'\xi$ for some $c>0$ (in the queueing model $\sigma\sigma'$ has an explicit formula in terms of the arrival variability $c_a^2$, $p$, $P$, $v$; only the ellipticity is used). The limiting diffusion is

$$X_t=x+\int_0^t g(X_s)\,ds+\sigma B_t,$$

with $B$ a standard $d$-dimensional Brownian motion; $P_t(x,\cdot)$ is the law of $X_t$ started at $x$, and $\mu$ is a probability measure with $\int P_t(x,\cdot)\,\mu(dx)=\mu$ for all $t\ge0$.

The Euler–Maruyama chain with step $\eta>0$ is

$$\tilde X^\eta_{k+1}=\tilde X^\eta_k+g(\tilde X^\eta_k)\,\eta+\sqrt\eta\,\sigma\xi_{k+1},\qquad k\in\mathbb N_0,$$

with independent standard Gaussian $\xi_k$; its one-step kernel is $\tilde P_\eta(x,\cdot)=\mathcal N(x+\eta g(x),\eta\sigma\sigma')$, and $\tilde\mu_\eta$ denotes an invariant probability measure of it. The Wasserstein-1 distance is taken in its dual form

$$d_W(\mu_1,\mu_2)=\sup_{h\in\mathrm{Lip}(1)}\Big|\int h\,d\mu_1-\int h\,d\mu_2\Big|,$$

where $\mathrm{Lip}(1)$ is the set of $h$ with $|h(x)-h(y)|\le|x-y|$.

Formalization targets

Goal: Theorem 1.2(i)

There is $\eta_0>0$ such that for every $\varsigma\in(0,1)$ there is $C_\varsigma>0$, depending on $\varsigma$ but not on $\eta$, such that for every $\eta\in(0,\eta_0]$ the chain has a unique invariant probability measure $\tilde\mu_\eta$ and

$$d_W(\mu,\tilde\mu_\eta)\le C_\varsigma\,\eta^{\frac{1-\varsigma}{2}}.$$

The constants are left unspecified, so the goal survives any improvement of them.

Milestones

In the order the proof uses them:

1. (A.4)–(A.7): a $C^2$ Lyapunov function $V$ with $\hat c_1|y|^2\le V(y)\le\hat C|y|^2+\hat C'$, $|\nabla V(y)|\le C(1+|y|)$ and $\mathcal AV\le-c_1V+\breve c_1$, where $\mathcal A$ is the generator.
2. Lemma A.1: moment bounds $\mathbb EV^\ell(X^x_t)\le e^{-c_1t}V^\ell(x)+\breve c_\ell(1-e^{-c_1t})/c_1$, and exponential ergodicity of the diffusion in $d_W$ and in weighted total variation.
3. Lemma A.2: $\tilde P_\eta f$ is Lipschitz for bounded measurable $f$ (strong Feller), and the chain is irreducible.
4. Lemma A.3: the one-step error $\mathbb E|X^x_\eta-\tilde X^{\eta,x}_1|^{2\ell}\le\tilde C_\ell(1+V^\ell(x))\eta^{3\ell}$.
5. Lemma A.4: the drift condition $\tilde P_\eta V^\ell\le(e^{-c_1\eta}+\tilde C_\ell\eta^{3/2})V^\ell+\check K_\ell$ for the chain.
6. Proposition 2.1: existence, uniqueness and exponential ergodicity of $\tilde\mu_\eta$, and $\tilde\mu_\eta(|\cdot|^\ell)\le C$ uniformly in $\eta$.
7. Lemma 2.2: for the solution $f$ of the Stein equation $\mathcal Af=h-\mu(h)$ with $h\in\mathrm{Lip}_0(1)$, the polynomial bounds on $f,\nabla f,\nabla^2f$ and the local $(1-\varsigma)$-Hölder bound on $\nabla^2 f$.
8. (2.7)–(2.8): the second-order Taylor remainder in Stein's argument is $O(\eta^{(3-\varsigma)/2})$.

Significance

The theorem gives a computable approximation of the stationary distribution of the diffusion limit, with an explicit rate. Combined with the convergence of the scaled queue's stationary distribution to $\mu$ at rate $n^{(\varsigma-1)/2}$ for $M/Ph/n+M$ queues (Braverman and Dai), it yields an approximation of the steady state of the queue itself, by simulating the chain for $O(n\log n)$ steps (Remark 2.3 of the paper). Proposition 2.1 also shows that the explicit Euler scheme does not blow up for this non-gradient, non-differentiable drift, which is not automatic: the Euler scheme of an ergodic SDE may blow up.

The result is proved in the paper; nothing here is open mathematically, though Lemma 2.2 is imported from Gurvich's regularity theory rather than proved there. To our knowledge none of it is formalized. A formalization provides the first machine-checked Wasserstein bound between the invariant measures of a multidimensional diffusion and its discretization, and on the way the Lyapunov and Harris-type ergodicity theory of Euler–Maruyama chains.

Difficulty

Standard error analyses of Euler schemes for invariant measures assume a gradient drift, strong dissipativity, or several bounded derivatives; none holds here. The drift is not differentiable on the hyperplane $\{e'x=0\}$ and is dissipative only in the Lyapunov sense of Dieker and Gao, so contraction-based couplings do not apply. The Stein argument instead needs second derivatives of the Stein solution $f$ that are Hölder continuous with polynomially growing constants; this is where the loss $\eta^{-\varsigma/2}$ comes from, and it rests on elliptic regularity for an operator with only Lipschitz coefficients. The ergodicity of the chain must also hold with constants uniform in $\eta$, which the classical Lyapunov criteria do not give directly.

Formalization scope

• The state space is EuclideanSpace ℝ (Fin d), so $|\cdot|$ is the Euclidean norm; matrices act by the usual product, and $\|A\|_{\mathrm{op}}=\sup_{|u|=1}|Au|$.
• $\sigma$ is a general uniformly elliptic matrix, not the paper's explicit covariance; the general case contains the paper's.
• Brownian motion is the published EthierKurtz.IsStandardBrownian (time ℝ≥0). The diffusion is a pathwise integral equation with continuous paths, so no Itô integral is needed. The transition function $P_t$ is characterized by: on every probability space, every solution from $x$ has law $P_t(x,\cdot)$ at time $t$. $\mu$ is a hypothesis (an invariant probability measure of $P_t$), never assumed smooth or with moments.
• The chain's kernel is the image of stdGaussian under $z\mapsto x+\eta g(x)+\sqrt\eta\,\sigma z$; $k\in\mathbb N_0$ starts at $0$.
• $d_W$, $\|\cdot\|_{\mathrm{TV}}$ and $\|\cdot\|_{\mathrm{TV},V}$ take values in $[0,\infty]$ and equal $+\infty$ unless the needed first moment (or $\int(1+V)$) is finite; so a bound on $d_W$ cannot hold because an integral of a non-integrable test function defaults to $0$. All moment bounds are lower Lebesgue integrals.
• "For any small enough $\varsigma$" is read as every $\varsigma\in(0,1)$; the step threshold $\eta_0$ is chosen before $\varsigma$, and $C_\varsigma$ before $\eta$.
• The Stein solution is $f=-\int_0^\infty(P_th-\mu(h))\,dt$ (the paper's (3.2)); Lemma 2.2 concludes that it is $C^2$ and solves $\mathcal Af=h-\mu(h)$, so its derivatives are never junk.
• Lemma A.2's gradient bound is stated in Lipschitz form, since $\tilde P_\eta f$ need not be differentiable where $g$ is not.

Trivializing encodings are excluded: the existence and uniqueness of $\tilde\mu_\eta$ is a conclusion, not a hypothesis, and the distance is infinite rather than zero when moments fail.

A complete development needs Brownian motion with Gaussian increments, a pathwise SDE solution theory for Lipschitz drift, Itô's formula (or Dynkin's formula) for $C^2$ functions of additive-noise diffusions, Harris-type ergodic theorems for Feller chains and diffusions, and elliptic regularity for the Stein equation. The Harris theorems and Wasserstein machinery are reusable well beyond this mission. Contributions to any milestone, and to general infrastructure for Euler–Maruyama chains, are welcome.

Selected references

• X. H. Jin, G. Pang, L. Xu, X. Xu, An approximation to the invariant measure of the limiting diffusion of G/Ph/n+GI queues in the Halfin–Whitt regime and related asymptotics, arXiv:2209.07361v1, 2022. https://arxiv.org/abs/2209.07361
• A. B. Dieker, X. Gao, Positive recurrence of piecewise Ornstein–Uhlenbeck processes and common quadratic Lyapunov functions, Annals of Applied Probability 23(4): 1291–1317, 2013.
• J. G. Dai, S. He, T. Tezcan, Many-server diffusion limits for G/Ph/n+GI queues, Annals of Applied Probability 20(5): 1854–1890, 2010.
• J. G. Dai, A. B. Dieker, X. Gao, Validity of heavy-traffic steady-state approximations in many-server queues with abandonment, Queueing Systems 78(1): 1–29, 2014.
• S. Browne, W. Whitt, Piecewise-linear diffusion processes, in J. Dshalalow (ed.), Advances in Queueing: Theory, Methods, and Open Problems, 463–480, 1995.
• I. Gurvich, Diffusion models and steady-state approximations for exponentially ergodic Markovian queues, Annals of Applied Probability 24(6): 2527–2559, 2014.
• A. Braverman, J. G. Dai, Stein's method for steady-state diffusion approximations of M/Ph/n+M systems, Annals of Applied Probability 27(1): 550–581, 2017.
• M. Hairer, J. C. Mattingly, Spectral gaps in Wasserstein distances and the 2D stochastic Navier–Stokes equations, Annals of Probability 36(6): 2050–2091, 2008.

Bibliographic data as in the reference list of arXiv:2209.07361v1, pp. 35–37.

Motivation

In many sequential decision problems the state of a system is not observed for free. A controller monitoring a machine, a patient or an inventory must pay to inspect it, and between inspections it acts on a belief rather than on the state. The observation cost model (OCM) makes this precise for a finite Markov chain: at each time the controller either pays a fixed cost $c_{\mathrm{obs}}>0$ to observe the state and choose a new action, or keeps the previous action without observing. Reisinger and Tam (arXiv:2201.07908v3, published in Mathematics of Operations Research, 2025) show that, after augmenting the state with the time elapsed since the last observation, the dynamic programming equation of the discounted OCM is a discrete quasi-variational inequality (QVI): an obstacle problem whose obstacle depends on the unknown itself.

Earlier work on the OCM (Anderson and Friedman, 1977; Hajek, Mitzel and Yang, 2008) treated fixed dynamics of the underlying chain. Discrete QVIs with interconnected obstacles and their penalty approximation were analysed by Reisinger and Zhang (2019) for switching problems, with the identity matrix in place of the transition operator $Q_n$. The paper recasts the OCM QVI in an abstract form and proves it is well-posed.

Setting

Let $X$ be a finite state space and $A$ a finite nonempty action set. The unknown is a function $u:\mathbb N\times X\to\mathbb R^A$ with components $u_a(n,x)=u^n_{a,x}$. A driver is a family $F_a:\mathbb N\times X\times\mathbb R^A\times\mathbb R^{\mathbb N\times X}\to\mathbb R$, $a\in A$; its third argument is the vector $u(n,x)$ and its fourth the whole function $u_a$. Given a sequence $\{Q_n\}$ of substochastic matrices on $X$ (nonnegative entries, row sums at most $1$) and a vector $c\in\mathbb R^{X\times A}$, the inspection operator is

$$\mathcal Mu(n,x)=\big(Q_n\,\overline{u^1-c}\big)_x,\qquad \big(\overline{u^1-c}\big)_y=\max_{a\in A}\big(u^1_{a,y}-c_{a,y}\big).$$

It reads $u$ only at time $1$ and couples all actions. The QVI (3.1) is

$$\min\big\{F_a(n,x,u(n,x),u_a),\;u_a(n,x)-\mathcal Mu(n,x)\big\}=0,\qquad n\in\mathbb N,\ x\in X,\ a\in A.$$

Assumption 3.2 asks for a constant $\beta>0$ with $F_a(n,x,r,u_a+\theta)-F_a(n,x,s,v_a)\ge\beta\theta$ whenever $u\le v$ are bounded and $\theta=r_a-s_a=\max_j(r_j-s_j)\ge0$, together with boundedness in $x$ and uniform continuity in $r$. Assumption 3.12 asks that $F_a(n,x,r,\phi)$ depend on $\phi$ only through finitely many values. For the OCM, $F_a(n,x,r,\phi)=r_a-\gamma\phi(n+1,x)-(P^n_ar_a)_x$.

For the computation the paper truncates time to $\{1,\dots,N\}$, giving a finite QVI (3.7) for a map $\tilde F:\mathbb R^{N\times L\times d}\to\mathbb R^{N\times L}$ with the monotonicity property (3.5), and replaces it by the penalised problem (3.8), $\tilde F_a(u)-\rho\,\pi(\mathcal Mu-u_a)=0$, with a continuous non-decreasing penalty $\pi$ vanishing exactly on $(-\infty,0]$.

Formalization targets

Goal: well-posedness (Theorem 3.13)

Under Assumptions 3.2 and 3.12, a uniform bound $\sup_{a,n,x}|F_a(n,x,0,0)|<\infty$, and $c_{a,x}>0$:

$$\exists!\,u:\mathbb N\times X\to\mathbb R^A\ \text{bounded with}\ \min\{F_a(n,x,u(n,x),u_a),\,u_a-\mathcal Mu\}=0\ \text{everywhere}.$$

Milestones

• Lemma 3.5: the restriction of $F$ to $\{1,\dots,N\}$ satisfies (3.5) with the same $\beta$.
• Proposition 3.10: for every $\rho\ge0$, a subsolution of (3.8) lies below a supersolution.
• Lemma 3.11: for $c\ge0$, every solution of (3.8) satisfies $\|u^\rho\|_\infty\le\|\tilde F(0)\|_\infty/\beta$.
• Theorem 3.13, truncated part: (3.8) has a unique solution $u^{\rho,N}$, which increases to a limit $u^N$ as $\rho\to\infty$; for $c>0$, $u^N$ solves (3.7).
• Proposition 3.8: comparison for the truncated QVI (3.7) when $c>0$.
• Proposition 3.3: comparison for bounded sub- and supersolutions of (3.1) when $c>0$.
• Corollary 3.4: bounded solutions of the OCM QVI (2.11) are unique.

Significance

The well-posedness theorem identifies the value function of the discounted observation cost model as the unique bounded solution of a QVI, which is what justifies computing it: truncating time and penalising the obstacle gives a sequence of finite nonlinear systems whose solutions converge to it, and these systems are solved in practice by semismooth Newton iterations. The abstract framework covers other models with a non-local, action-coupling obstacle, such as discretised switching problems.

The results are proved in the paper, partly by reference: existence of the penalised solutions and the convergence $u^{\rho,N}\to u^N$ are delegated to Reisinger and Zhang (2019, Theorems 2.5 and 2.6), and Lemma 3.11 is stated as "a direct adaptation" of a lemma there. None of the results is formalized. A formal development makes the delegated steps explicit and fixes the hypotheses under which the statements hold (see Formalization scope).

Difficulty

The comparison principle (Proposition 3.3) is a maximum principle on an infinite index set: the supremum of $u-v$ need not be attained, so the argument works along maximising sequences and needs the uniform continuity of the driver in $r$. In the obstacle branch, the inspection operator transports the maximum to time $1$, and excluding that case uses $c>0$ through a chain argument over the finite set $A\times X$; with $c=0$ comparison fails.

Existence is harder. A direct fixed-point argument fails because the driver is only monotone, not a contraction, and the obstacle depends on the unknown. The paper's route goes through two limits, $\rho\to\infty$ and $N\to\infty$, each requiring bounds uniform in the other parameter and a passage to the limit in a nonlinear, non-local equation; the $N\to\infty$ limit uses Assumption 3.12 to localise the dependence on $\phi$.

Formalization scope

Everything is finite-dimensional or a function on the discrete set $\mathbb N\times X$; no measure theory is involved. Unknowns are curried as u n x a $=u_a(n,x)$, costs as c a x. Time in (3.1) ranges over all of $\mathbb N$, including $0$, and $u^1$ is u 1. The truncated time set $\{1,\dots,N\}$ is Fin N with the index $k$ standing for time $k+1$, so $u^1$ is u 0. Substochastic means nonnegative entries with row sums at most $1$; $c>0$ and $c\ge0$ are componentwise; norms are sup norms, written out entrywise. Bounded sets of $r$ in Assumption 3.2.2 are cubes.

Two additions to the paper's hypotheses are made, each because the statement is false without it. The goal assumes $\sup_{a,n,x}|F_a(n,x,0,0)|<\infty$: Assumption 3.2.2 only bounds $F$ for each fixed $n$, and the driver $r_a-\gamma\phi(n+1,x)-n$ satisfies both assumptions yet has no bounded solution. Corollary 3.4 is stated for bounded solutions of (2.11). The truncated QVI is stated without Problem 3.6's time boundary condition $u_a(N,\cdot)=\mathcal Mu(N,\cdot)$, which (3.8) does not encode. The continuity of the restricted driver claimed in Lemma 3.5 is not part of that milestone; the truncated part of Theorem 3.13 takes continuity of an abstract $\tilde F$ as a hypothesis.

The $N\to\infty$ clause of Theorem 3.13 ("$u=\lim_{N\to\infty}u^N$, monotone from below") is not stated: the paper does not define the truncated operator beyond time $N$, and with the zero extension the monotonicity in $N$ fails. The goal keeps the clause's well-defined content, existence and uniqueness of the bounded solution. A goal asserting only uniqueness (Proposition 3.3), or existence with the uniform bound hidden inside a definition, would be a different theorem; the bound is an explicit hypothesis.

A complete development needs solvability of finite-dimensional monotone nonlinear systems (the existence result cited from Reisinger and Zhang), monotone convergence in finite dimensions, and maximum-principle arguments on $\mathbb N\times X$. The comparison principles and the a-priori bound are reusable for other monotone schemes. Proofs of any milestone, and alternative existence arguments for the goal, are welcome.

Selected references

• C. Reisinger and J. Tam, Markov decision processes with observation costs: framework and computation with a penalty scheme, arXiv:2201.07908v3, 2023; Mathematics of Operations Research, 2025. https://arxiv.org/abs/2201.07908v3
• C. Reisinger and Y. Zhang, A penalty scheme for monotone systems with interconnected obstacles: convergence and error estimates, SIAM Journal on Numerical Analysis 57(4):1625–1648, 2019 (reference [28] of the paper).
• R. F. Anderson and A. Friedman, Optimal inspections in a stochastic control problem with costly observations, Mathematics of Operations Research 2(2):155–190, 1977. https://doi.org/10.1287/moor.2.2.155
• B. Hajek, K. Mitzel and S. Yang, Paging and registration in cellular networks: jointly optimal policies and an iterative algorithm, IEEE Transactions on Information Theory 54:608–622, 2008.

Why the smoothness of an image function matters

The alternating direction method of multipliers (ADMM) solves problems of the form $\min\{f(x)+g(z)\mid Ax+Bz=b\}$ and is a standard tool in signal processing, statistics and distributed optimization (Boyd et al., 2011). For convex problems ADMM is known to be the Douglas–Rachford splitting (DRS) applied to a dual problem. Themelis and Patrinos (arXiv:1709.05747v4, SIAM J. Optim. 30(1), 2020) show a primal equivalence that holds without convexity: one ADMM step is one DRS step on the problem

$$\min_{s\in\mathbb R^p}\ (Af)(s)+(Bg)(b-s),$$

where $(Af)$ and $(Bg)$ are image functions (their Theorem 5.5). Their nonconvex DRS theory (Theorem 4.3) requires the first summand $\varphi_1=(Af)$ to have a Lipschitz continuous gradient. Convergence guarantees for nonconvex ADMM therefore rest on knowing when $(Af)$ is smooth and with which constants. Section 5.4.2 of the paper answers this question in Theorem 5.13, the target of this mission.

Setting

Let $A\in\mathbb R^{p\times n}$ be a matrix, viewed as a linear map $\mathbb R^n\to\mathbb R^p$ with transpose $A^\top$ and spectral norm $\|A\|$, and let $f:\mathbb R^n\to\mathbb R$. The image function (Definition 5.1) is

$$(Af)(s)=\inf\{f(x)\mid Ax=s\},\qquad s\in\mathbb R^p,$$

the least value of $f$ on the fibre $A^{-1}(s)$. It can be $+\infty$ (empty fibre) or $-\infty$ (unbounded below). Write $X(s)=\operatorname*{arg\,min}\{f(x)\mid Ax=s\}$ for the set of minimizers on the fibre.

A differentiable $h$ is $L$-smooth if $\|\nabla h(x)-\nabla h(y)\|\le L\|x-y\|$, and $\sigma$-hypoconvex if $h-\frac\sigma2\|\cdot\|^2$ is convex; $C^{1,1}(\mathbb R^n)$ is the class of $L$-smooth functions (§2.2). For $\sigma>0$ hypoconvexity is $\sigma$-strong convexity; for $\sigma<0$ it allows nonconvexity.

The paper introduces smoothness relative to a matrix (Definition 5.12). A differentiable $f$ belongs to $C^{1,1}_A(\mathbb R^n)$ with constants $|\sigma_{f,A}|\le L_{f,A}$ if

$$\sigma_{f,A}\|A(x-y)\|^2\le\langle\nabla f(x)-\nabla f(y),x-y\rangle\le L_{f,A}\|A(x-y)\|^2$$

whenever $\nabla f(x),\nabla f(y)\in\operatorname{range}A^\top$. The condition only concerns points whose gradient is orthogonal to $\ker A$, which includes every minimizer of $f$ on a fibre.

Finally, $\sigma_+(M)$ denotes the smallest nonzero singular value of a nonzero matrix $M$ (p. 5).

Formalization targets

Goal: Theorem 5.13

Let $A$ be surjective and $f$ lower semicontinuous, and suppose that for some $\beta\ge0$ the function $f+\frac\beta2\|A\cdot-s\|^2$ is level bounded for every $s$. Then $(Af)$ is real-valued, $L_{(Af)}$-smooth and $\sigma_{(Af)}$-hypoconvex on $\mathbb R^p$ in each of three cases:

1. $f\in C^{1,1}_A(\mathbb R^n)$: $\quad L_{(Af)}=L_{f,A}$, $\ \sigma_{(Af)}=\sigma_{f,A}$;
2. $f\in C^{1,1}(\mathbb R^n)$ and $X(s)$ is single valued and $M$-Lipschitz:

$$L_{(Af)}=L_fM^2,\qquad \sigma_{(Af)}=\begin{cases}\sigma_f/\|A\|^2&\sigma_f\ge0,\\ \sigma_fM^2&\sigma_f<0;\end{cases}$$

3. $f\in C^{1,1}(\mathbb R^n)$ is convex: $\quad L_{(Af)}=L_f/\sigma_+(A^\top A)$, $\ \sigma_{(Af)}=\sigma_f/\|A\|^2$.

Milestones

In the order the proof uses them: properness of image functions (Proposition 5.2(i)); lower semicontinuity of $(Af)$ and nonemptiness of $X(s)$; strict continuity of $(Af)$ and nonemptiness of its limiting subdifferential; the subgradient transfer $A^\top\hat\partial(Af)(\bar s)\subseteq\hat\partial f(\bar x)$ (Proposition 5.3); the identity (5.7) relating subgradients of $(Af)$ to gradients of $f$ at fibre minimizers; strong convexity of image functions (Proposition 5.4); and the subdifferential characterization of smoothness (Lemma 2.1).

Significance

Theorem 5.13 is what makes the DRS convergence results of the paper apply to ADMM in practice: Theorem 5.6 of the paper assumes $(Af)$ is smooth, and this theorem gives conditions on $f$ and $A$ under which that assumption holds, with explicit constants that enter the admissible penalty parameters. Case (i) is weaker than requiring $f$ itself to be $C^{1,1}$: the paper's example $f(x,y)=\frac12x^2y^2$, $A=[1\ 0]$ satisfies it with $L_{f,A}=0$ although $\nabla f$ is not globally Lipschitz. Case (iii) recovers the convex setting with a constant governed by the conditioning of $A$. Propositions 5.3 and 5.4 and Lemma 2.1 are general tools of variational analysis that apply beyond this theorem.

The result is proved in the paper, with parts of the argument delegated to Rockafellar–Wets (Thms. 1.32, 9.13, 10.13) and to Nesterov's cocoercivity theorem. No machine-checked version is known to exist. A formalization would also need image functions, regular and limiting subdifferentials on $\mathbb R^p$, and parametric minimization, none of which Mathlib currently provides.

Difficulty

The image function is defined by an infimum, so smoothness of $(Af)$ cannot be read off a formula. The proof goes through the limiting subdifferential: it shows that $(Af)$ is lower semicontinuous and strictly continuous, so that $\partial(Af)$ is nonempty, relates each subgradient $y$ of $(Af)$ to a gradient $\nabla f(x)$ at a fibre minimizer, and then applies a subdifferential characterization of $C^{1,1}$ functions (Lemma 2.1). The obvious shortcut, an envelope or implicit-function argument that differentiates through the minimizer map $X(s)$, fails: $X(s)$ need not be single valued or continuous in cases (i) and (iii). Limiting subgradients also require a compactness argument (uniform level boundedness) to pass from regular subgradients at nearby points to minimizers on the limit fibre.

Formalization scope

Spaces are EuclideanSpace ℝ (Fin n), matrices are continuous linear maps with ContinuousLinearMap.adjoint as transpose and the operator norm as $\|A\|$, and "surjective" is Function.Surjective. $f$ is real-valued; the image function imageFn is an EReal infimum over the fibre, and the conclusion asserts that it equals a real function $F$ that is $L$-smooth and $\sigma$-hypoconvex. Regular and limiting subdifferentials are the published NonconvexSplitting.Shared.IsRegularSubgrad and LimitingSubdiff. The three cases are three separate implications, each with the constants printed in the theorem. In case (iii), $\sigma_+(A^\top A)$ is passed as a number $c$ with the hypothesis that $c$ is the least value of $\|A^\top y\|^2$ over unit vectors $y$, which equals the smallest nonzero eigenvalue of $A^\top A$ for surjective $A$. For $p=0$ this set is empty and case (iii) is vacuous.

The goal does not assume that $(Af)$ is real-valued, differentiable or strictly continuous, nor that fibre minimizers exist outside case (ii). A statement making any of these assumptions would be trivial or circular. Definition 5.12 keeps its range guard; dropping the guard would turn case (i) into ordinary $C^{1,1}$ smoothness.

Three points differ from the printed text, and each is disclosed in the item concerned:

• Lemma 2.1 is stated for lower semicontinuous $h$. On $\mathbb R$, the step function $h=-1$ on $(-\infty,0)$ and $0$ on $[0,\infty)$ has $\partial h\equiv\{0\}$ but is not differentiable. Theorem 5.13 applies the lemma to the lsc function $(Af)$.
• Strict continuity of $(Af)$ is stated for $C^1$ functions $f$. The page's "since $f$ is differentiable" is not enough: $f(x)=x^2\sin(1/x^2)+2x^2$ with $A=1$ is a counterexample.
• (5.7) is stated for limiting subgradients, which is what the proof needs.

The definitions of this mission (image function, relative smoothness, strong convexity of extended-valued functions) are reusable for any analysis of ADMM or of marginal functions. Proofs of the milestones, of Lemma 2.1 as a standalone result, and supporting variational-analysis lemmas are all welcome.

Selected references

• A. Themelis, P. Patrinos, Douglas–Rachford splitting and ADMM for nonconvex optimization: tight convergence results, SIAM J. Optim. 30(1):149–181, 2020. arXiv:1709.05747v4, doi:10.1137/18M1163993
• R. T. Rockafellar, R. J.-B. Wets, Variational Analysis, Grundlehren 317, Springer, 1998. doi:10.1007/978-3-642-02431-3
• Y. Nesterov, Lectures on Convex Optimization, Springer, 2018. doi:10.1007/978-3-319-91578-4
• H. H. Bauschke, P. L. Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Spaces, 2nd ed., Springer, 2017. doi:10.1007/978-3-319-48311-5
• S. Boyd, N. Parikh, E. Chu, B. Peleato, J. Eckstein, Distributed optimization and statistical learning via the alternating direction method of multipliers, Found. Trends Mach. Learn. 3(1):1–122, 2011. doi:10.1561/2200000016

Motivation

The expected error of the randomized SVD with a $k\times(k+p)$ Gaussian test matrix is $\bigl(1+\tfrac{k}{p-1}\bigr)^{1/2}$ times the optimal error (Halko–Martinsson–Tropp 2011, Theorem 10.5). The factor $\tfrac{k}{p-1}$ is $\mathbb E\|G^{\dagger}\|_F^2$ for a $k\times(k+p)$ standard Gaussian matrix $G$: the quality of every Gaussian sketching method is governed by inverse moments of a Gaussian matrix — moments of its pseudoinverse, equivalently of the inverse of the Wishart matrix $GG^{\mathsf T}$. The expected spectral error uses $\mathbb E\|G^{\dagger}\|$ instead (Tropp–Webber 2023, Lemmas B.2–B.3), and deviation bounds use tail estimates for both norms (HMT Proposition 10.4). None of these inverse moments is available in Mathlib. This mission proves them.

Setting

A standard Gaussian matrix $G\in\mathbb R^{r\times k}$ has independent $\mathcal N(0,1)$ entries; its law is $\gamma_{r,k}$. Throughout, $G$ is wide or square, $r\le k$, so that $GG^{\mathsf T}\in\mathbb R^{r\times r}$ is almost surely invertible; $GG^{\mathsf T}$ is a Wishart matrix with $k$ degrees of freedom and $(GG^{\mathsf T})^{-1}$ is inverse-Wishart.

The pseudoinverse of a full-row-rank matrix is $G^{\dagger}=G^{\mathsf T}(GG^{\mathsf T})^{-1}$. Its norms are determined by the inverse Wishart matrix:

$$\|G^{\dagger}\|_F^2=\operatorname{tr}(GG^{\mathsf T})^{-1},\qquad \|G^{\dagger}\|^2=\|(GG^{\mathsf T})^{-1}\|=\sigma_{\min}(G)^{-2},$$

where $\|\cdot\|_F$ is the Frobenius norm, $\|\cdot\|$ the spectral norm (largest singular value) and $\sigma_{\min}(G)$ the smallest of the $r$ singular values of $G$. The letter $e$ denotes Euler's number.

Formalization targets

Goal — Tropp–Webber inverse spectral moment (Lemma B.3)

For integer $1\le p\le18$ and $k\ge r+2p$,

$$\Bigl(\mathbb E\,\bigl\|(GG^{\mathsf T})^{-1}\bigr\|^p\Bigr)^{1/p}\le\frac{e^2(k+r)}{2(k-r)^2}.$$

The case $p=1$ is $\mathbb E\|G^{\dagger}\|^2\le e^2(k+r)/(2(k-r)^2)$, the spectral inverse moment at the heart of expected spectral-norm error bounds. It is the goal because it is the one inverse moment that is not an exact formula: its proof passes through Edelman's bound on the density of $\lambda_{\min}(GG^{\mathsf T})$ and an explicit integration, and it implies the expectation bound of HMT Proposition 10.2 by Jensen's inequality.

Milestones

• Inverse-Wishart mean. $\mathbb E(GG^{\mathsf T})^{-1}=\frac{1}{k-r-1}I_r$ for $r\le k-2$ (Tropp–Webber Lemma B.2).
• Frobenius inverse moment. $\mathbb E\|G^{\dagger}\|_F^2=\frac{r}{k-r-1}$ for $k-r\ge2$ (HMT Proposition 10.2 (10.3) and A.5).
• Second-order inverse-Wishart moments. $\mathbb E\|G^{\dagger}\|_F^4=\frac{r^2(k-r)-2r(r-1)}{(k-r)(k-r-1)(k-r-3)}$ and $\mathbb E\|(GG^{\mathsf T})^{-1}\|_F^2=\frac{r(k-1)}{(k-r)(k-r-1)(k-r-3)}$ for $r\le k-4$ (Tropp–Webber Lemma B.2).
• Chen–Dongarra tail. $\mathbb P\{\|G^{\dagger}\|>t\}\le\frac{1}{\sqrt{2\pi(k-r+1)}}\bigl(\frac{e\sqrt k}{k-r+1}\bigr)^{k-r+1}t^{-(k-r+1)}$ for $k\ge r\ge2$, $t>0$ (HMT Proposition A.3).
• Expected spectral norm of the pseudoinverse. $\mathbb E\|G^{\dagger}\|\le\frac{e\sqrt k}{k-r}$ for $r\ge2$, $k\ge r+1$ (HMT Proposition 10.2 (10.4) and A.4).
• Frobenius tail. $\mathbb P\{\|G^{\dagger}\|_F^2>\frac{12r}{k-r}t\}\le4t^{-(k-r)/2}$ for $k-r\ge4$, $t\ge1$ (HMT Theorem A.6).

Significance

The results. Combined with the second-moment identity and block independence of Mission I, the Frobenius inverse moment gives the expected Frobenius error of the randomized SVD and of generalized Nyström; the spectral inverse moment gives the expected spectral error and the core-noise bounds of truncated sketching methods; the two tail bounds give the large-deviation statements of HMT §10.3. The second-order inverse-Wishart moments quantify the variance of $\|G^{\dagger}\|_F^2$, which Tropp–Webber use for concentration of the randomized SVD error.

Formalizing them. The inverse-Wishart mean and the second-order moments are classical multivariate-statistics formulas (Muirhead; von Rosen) with no Mathlib counterpart. The Chen–Dongarra tail requires the exact distribution of $\sigma_{\min}$, or at least Edelman's density bound; the Frobenius tail is HMT's own argument through the bidiagonal (Householder) model of a Gaussian matrix and chi-square moment estimates. Every statement here is a proved theorem in the literature; none has a machine-checked proof.

Difficulty

There is no Wishart distribution in Mathlib. The direct approach — compute the density of $GG^{\mathsf T}$ and integrate — needs the Wishart density, a change of variables on the cone of positive-definite matrices, and the matrix-variate gamma function. The classical route to the mean avoids this: by rotational invariance $(GG^{\mathsf T})^{-1}$ has the same law as $Q(GG^{\mathsf T})^{-1}Q^{\mathsf T}$ for every orthogonal $Q$, so its mean is a multiple of the identity, and the trace is $\sum_i 1/\chi^2_{k-r+1}$-type quantities whose means are elementary; making "the diagonal entries of $(GG^{\mathsf T})^{-1}$ are reciprocals of $\chi^2_{k-r+1}$ variables" rigorous requires the Gram–Schmidt (QR) decomposition of a Gaussian matrix and the independence of its factors. The spectral inverse moment needs a lower-tail bound for $\lambda_{\min}(GG^{\mathsf T})$ of the exact order $t^{(k-r+1)/2}$; a union bound over rows loses a factor $r$ and does not give the stated constant. The naive attempt to bound $\mathbb E\|G^{\dagger}\|^2$ by $\mathbb E\|G^{\dagger}\|_F^2$ loses a factor $r$ as well.

Formalization scope

$G$ is sampled from the product measure $\gamma_{r,k}$ on arrays; $(\cdot)^{-1}$ is Mathlib's total matrix inverse, equal to $0$ on the null set where $GG^{\mathsf T}$ is singular, and $G^{\dagger}$ is defined as $G^{\mathsf T}(GG^{\mathsf T})^{-1}$ (pinvR). Every inequality is stated together with integrability of its integrand, so the Lean convention that a non-integrable function integrates to $0$ cannot make a bound vacuous; identities with nonzero right-hand sides force integrability on their own but are also stated with it. The inverse-Wishart mean is stated entrywise. Tail bounds are inequalities between a measure in $[0,\infty]$ and $\operatorname{ofReal}$ of the real bound. Real powers with non-integer exponents are Mathlib's rpow; all bases are positive under the hypotheses. The spectral norm is Mathlib's $\ell_2$ operator norm. The sizes follow Tropp–Webber ($G\in\mathbb R^{r\times k}$); HMT's $k\times(k+p)$ and $m\times n$ conventions are translated in each item's source field.

Welcome contributions: a Wishart/inverse-Wishart library (QR decomposition of a Gaussian matrix, Bartlett decomposition, chi-square moments $\mathbb E\,\Xi^{-q}=\Gamma(k/2-q)/(2^q\Gamma(k/2))$), and the distribution of $\sigma_{\min}$ of a Gaussian matrix (Edelman).

Selected references

• N. Halko, P.-G. Martinsson, J. A. Tropp, Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions, SIAM Review 53(2), 217–288, 2011. https://arxiv.org/abs/0909.4061
• J. A. Tropp, R. J. Webber, Randomized algorithms for low-rank matrix approximation: design, analysis, and applications, 2023. https://arxiv.org/abs/2306.12418
• R. Vershynin, Introduction to the non-asymptotic analysis of random matrices, in Compressed Sensing: Theory and Applications, Cambridge University Press, 2012. https://arxiv.org/abs/1011.3027
• L. Isserlis, On a formula for the product-moment coefficient of any order of a normal frequency distribution in any number of variables, Biometrika 12, 134–139, 1918. https://doi.org/10.1093/biomet/12.1-2.134
• S. Dasgupta, A. Gupta, An elementary proof of a theorem of Johnson and Lindenstrauss, Random Structures & Algorithms 22, 60–65, 2003. https://doi.org/10.1002/rsa.10073
• M. Rudelson, R. Vershynin, The Littlewood–Offord problem and invertibility of random matrices, Advances in Mathematics 218, 600–633, 2008. https://arxiv.org/abs/math/0703503
• Y. Chen, J. Dongarra, Condition numbers of Gaussian random matrices, SIAM J. Matrix Anal. Appl. 27(3), 603–620, 2005. https://doi.org/10.1137/040616413
• A. Edelman, Eigenvalues and condition numbers of random matrices, SIAM J. Matrix Anal. Appl. 9(4), 543–560, 1988. https://doi.org/10.1137/0609045
• Y. Gordon, Some inequalities for Gaussian processes and applications, Israel J. Math. 50, 265–289, 1985. https://doi.org/10.1007/BF02759761
• K. R. Davidson, S. J. Szarek, Local operator theory, random matrices and Banach spaces, Handbook of the Geometry of Banach Spaces, Vol. 1, 317–366, 2001. https://doi.org/10.1016/S1874-5849(01)80010-3
• M. Ledoux, The Concentration of Measure Phenomenon, AMS Mathematical Surveys and Monographs 89, 2001. https://doi.org/10.1090/surv/089

Shared definitions: the items of this mission import Definitions.Def_GaussianMatrix_basic, a draft of Mission I; submit that proposal first so the import resolves at compile time.

Series: I. Rotation Invariance and Moment Identities · II. Inverse Moments of the Pseudoinverse · III. Extreme Singular Values and Concentration. All three missions share the definition bundle GaussianMatrix_basic, published with Mission I, and the namespace GaussianMatrix.

Motivation

The alternating direction method of multipliers (ADMM) and the closely related augmented Lagrangian method are workhorses for linearly constrained optimization: they split a problem into cheap primal steps and a dual (multiplier) update, and they scale to large and distributed instances (Boyd et al., 2011). For convex objectives their convergence theory is classical (Eckstein–Bertsekas, 1992). Many applications, however, have nonconvex objectives: matrix factorization, sparse principal component analysis, and problems from signal processing and machine learning. There the simple inexact augmented Lagrangian method, which takes one projected-gradient step on the augmented Lagrangian before each multiplier update, can oscillate or diverge, and an example of this is given in §6 of the paper.

Zhang and Luo (SIAM J. Optim. 30(3), 2020) stabilize the method with a proximal term centred at an exponentially averaged copy of the primal iterates. They prove global convergence to stationary points for a smooth nonconvex objective, under linear equality constraints and box constraints, and an $O(1/\varepsilon^2)$ iteration complexity, without assuming that $A$ has full row rank. Earlier analyses of nonconvex ADMM, for example Hong–Luo–Razaviyayn, 2016 and Boţ–Nguyen, 2020, rely on other structure, such as a smooth last block or Kurdyka–Łojasiewicz properties. This mission formalizes the convergence theorem.

Setting

Let $f:\mathbb R^n\to\mathbb R$ be differentiable, $A\in\mathbb R^{m\times n}$, $b\in\mathbb R^m$, and let $P=\{x\in\mathbb R^n:\ \ell_i\le x_i\le u_i\}$ be a box with $\ell_i<u_i$. The problem is

$$\min_x f(x)\quad\text{subject to}\quad Ax=b,\ x\in P. \tag{1.1}$$

A point $x$ is stationary if there are multipliers $y\in\mathbb R^m$ and $\mu,\nu\in\mathbb R^n$ with

$$\nabla f(x)+A^\top y-\mu+\nu=0,\quad Ax=b,\quad x\in P,\quad \mu,\nu\ge 0,\quad \mu_i(\ell_i-x_i)=0,\quad \nu_i(x_i-u_i)=0 .$$

$X^*$ is the set of stationary points and $W$ the set of primal-dual pairs $(x,y)$ for which such $\mu,\nu$ exist. Strict complementarity means that for every such tuple and every $i$, exactly one of $\mu_i$ and $\ell_i-x_i$ is zero and exactly one of $\nu_i$ and $x_i-u_i$ is zero.

Assumption 2.2 asks for three things: (a) $0$ lies in the relative interior of $\{Ax-b: x\in P\}$; (b) strict complementarity holds; (c) $\|\nabla f(x)-\nabla f(x')\|\le L\|x-x'\|$ on $P$ for some $L>0$. A constant $\gamma$, possibly negative, satisfies $\langle\nabla f(x)-\nabla f(x'),x-x'\rangle\ge\gamma\|x-x'\|^2$ on $P$.

With constants $\Gamma>0$ and $p>0$, the proximal augmented Lagrangian is

$$K(x,z;y)=f(x)+y^\top(Ax-b)+\frac\Gamma2\|Ax-b\|^2+\frac p2\|x-z\|^2 .$$

Algorithm 2.2 starts from $x^0,z^0\in P$ and $y^0\in\mathbb R^m$, and for $t=0,1,\dots$ sets

$$y^{t+1}=y^t+\alpha(Ax^t-b),\qquad x^{t+1}=[x^t-c\nabla_xK(x^t,z^t;y^{t+1})]_+,\qquad z^{t+1}=z^t+\beta(x^{t+1}-z^t),$$

where $[\cdot]_+$ is the projection onto $P$, $\alpha>0$ and $0<\beta\le1$ are the dual and smoothing stepsizes, and $0<c<1/L_K$ with $L_K=L+p+\Gamma\|A\|^2$.

Formalization targets

Goal: Theorem 2.4

Assume Assumption 2.2, $p>-\gamma$, $p>0$, $\Gamma>0$ and $0<c<1/L_K$. There is $\alpha_0>0$ such that for every $\alpha\in(0,\alpha_0)$ and every $y^0$ there is $\beta_0>0$ such that for every $\beta\in(0,\min\{1,\beta_0\}]$, every run of Algorithm 2.2 from $y^0$ satisfies: $\{y^t\}$ is bounded,

$$\|x^{t+1}-x^t\|\to0,\qquad \operatorname{dist}(x^t,X^*)\to0,\qquad\operatorname{dist}(z^t,X^*)\to0,\qquad \operatorname{dist}((x^t,y^t),W)\to0,$$

and every limit point of $\{(x^t,y^t)\}$ lies in $W$. The thresholds are existential, as in the paper, which proves only that "sufficiently small" stepsizes work.

Milestones

The milestones follow the paper's proof:

1. weak duality (2.10);
2. the primal, dual and proximal descent lemmas (Lemmas 3.1–3.3);
3. the Lipschitz bounds of Lemma 3.5 and the primal error bounds (3.10)–(3.11);
4. the lower bound (3.5) on the potential $\phi^t=K(x^t,z^t;y^t)-2d(y^t,z^t)+2M(z^t)$;
5. dual boundedness (Lemma 3.4);
6. Lemma 3.6 and Corollaries 3.7–3.8;
7. the Hoffman bound (Proposition 3.9) and the local dual error bound (3.15);
8. the sufficient decrease

$$\phi^t-\phi^{t+1}\ge\frac1{8c}\|x^t-x^{t+1}\|^2+\frac\alpha2\|Ax(y^{t+1},z^t)-b\|^2+\frac p{6\beta}\|z^t-z^{t+1}\|^2 .\tag{3.37}$$

Here $d(y,z)$ and $x(y,z)$ are the minimum value and the minimizer of $K(\cdot,z;y)$ over $P$, and $M(z)$ and $x^*(z)$ are those of $f+\frac p2\|\cdot-z\|^2$ over the feasible set.

Companion: Theorem 4.2

Under the same hypotheses, for every run there is $B>0$ such that for every $t>0$ some $(x^{s+1},y^{s+1})$ with $s<t$ is a $B/\sqrt t$-solution: $\|Ax-b\|\le\varepsilon$ and some $v\in\nabla f(x)+A^\top y+\partial\iota_P(x)$ has $\|v\|\le\varepsilon$.

Significance

Theorem 2.4 is a global convergence guarantee for a single-loop, gradient-based primal-dual method on nonconvex problems with general linear constraints. Each iteration costs one gradient evaluation and one projection onto a box. Its hypotheses do not include full row rank of $A$ or a Kurdyka–Łojasiewicz property, and strict complementarity is generic by Proposition 2.3 of the paper. Theorem 4.2 turns the same argument into an $O(1/\varepsilon^2)$ iteration bound, which the paper reports as an improvement over the previously best known $O(1/\varepsilon^3)$.

The result is proved on paper and has not been machine-checked. A formalization would check the three-case argument for (3.37), whose printed Case 1 leaves a step implicit (it needs $\operatorname{dist}(z^t,X^*)\le\Delta$ but derives only a bound on $\|x^{t+1}-z^t\|$). It would also produce reusable components: a projected-gradient descent and error-bound library on boxes, Danskin's formula for the value of a strongly convex parametric problem, and Hoffman's bound for mixed linear systems.

Difficulty

None of the three terms $K$, $-d$ and $M$ decreases monotonically along the iteration, and the dual step can increase $K$ by $\alpha\|Ax^t-b\|^2$. A direct Lyapunov argument on $K$ therefore fails. The paper's potential cancels the dual increase against the dual ascent of $d$, but leaves a negative term $p\lambda\|x(y^{t+1},z^t)-x^*(z^t)\|^2$. Controlling that term needs a dual error bound, which holds only locally (near $X^*$ and for small dual residual) and only under strict complementarity. Away from that region, compactness of $P$ and the smallness of $\beta$ are used instead. The constants involved ($\Delta$, $\sigma_5$, $\delta(\cdot)$ from Corollaries 3.7–3.8) are not explicit, which is why the threshold for $\beta$ is not explicit either. Hoffman's bound and Danskin's theorem, both used in the proof, are not in Mathlib.

Formalization scope

Vectors live in EuclideanSpace ℝ (Fin n). $A$ is a continuous linear map, $A^\top$ is its adjoint, and the paper's $\sigma$ (largest singular value) is the operator norm $\|A\|$. The box projection is defined coordinatewise by the paper's formula (2.1). Assumption 2.2 is split into its three parts; the relative interior is Mathlib's intrinsicInterior, and "exactly one is zero" is an exclusive or. The Lipschitz bound and the constant $\gamma$ hold on $P$ only, as on the page.

The minimizers $x(y,z)$ and $x^*(z)$ are not built by choice. Statements take functions with the minimizer property, which exist because $P$ is compact and are unique when $p>-\gamma$. The values $d$, $M$ and the potential are computed at them. Distances are Metric.infDist, and on $\mathbb R^n\times\mathbb R^m$ this uses the maximum of the two Euclidean distances, an equivalent metric. A limit point is a cluster point of the sequence.

Standing assumptions of the paper that the Lean states explicitly: $\Gamma>0$, $p>0$, $p>-\gamma$, $0<c<1/L_K$, $0<\beta\le1$, $\alpha>0$. Additions relative to the page, each disclosed in the item: $x\in P$ in (2.10), $z\in P$ in Lemma 3.4 and in (3.15), nonemptiness of the polyhedron in Proposition 3.9, and the hypotheses of Theorem 2.4 for Theorem 4.2. The third constant of Lemma 3.5 is $\sigma/(p+\gamma)$, as in (3.12), not the printed $(p+\gamma)/\sigma$. Constants that the paper writes with a division by $\|A\|$ are multiplied out.

A trivializing formalization is ruled out. The stepsizes are quantified as "there is a threshold, then every smaller stepsize and every run". All six conclusions of Theorem 2.4 are stated. Its hypotheses are exactly the paper's. It does not mention the potential or the auxiliary minimizers. The sets $X^*$, $W$ and $Y^*(z)$ are nonempty under the hypotheses, so no distance to an empty set decides a statement. A sorry-free check shows that the hypotheses are jointly satisfiable, including strict complementarity, and that runs and minimizer selections exist.

Contributions welcome: proofs of the milestones in any order, general-purpose lemmas (projected-gradient error bounds, Danskin's theorem, Hoffman's bound) as separate theorems, and a repair of Case 1 of (3.37) if one is needed.

Selected references

• J. Zhang and Z.-Q. Luo, A proximal alternating direction method of multiplier for linearly constrained nonconvex minimization, SIAM J. Optim. 30(3):2272–2302, 2020. https://doi.org/10.1137/19M1242276
• S. Boyd, N. Parikh, E. Chu, B. Peleato and J. Eckstein, Distributed optimization and statistical learning via the alternating direction method of multipliers, Found. Trends Mach. Learn. 3(1):1–122, 2011. https://doi.org/10.1561/2200000016
• J. Eckstein and D. P. Bertsekas, On the Douglas–Rachford splitting method and the proximal point algorithm for maximal monotone operators, Math. Program. 55:293–318, 1992. https://doi.org/10.1007/BF01581204
• M. Hong, Z.-Q. Luo and M. Razaviyayn, Convergence analysis of alternating direction method of multipliers for a family of nonconvex problems, SIAM J. Optim. 26(1):337–364, 2016. https://doi.org/10.1137/140990309
• R. I. Boţ and D.-K. Nguyen, The proximal alternating direction method of multipliers in the nonconvex setting: convergence analysis and rates, Math. Oper. Res. 45(2):682–712, 2020. https://doi.org/10.1287/moor.2019.1008
• F. Facchinei and J.-S. Pang, Finite-Dimensional Variational Inequalities and Complementarity Problems, Springer, 2003. https://doi.org/10.1007/b97543

Motivation

A variational inequality asks for a point of a convex set at which a vector field points "inward" against every feasible direction. It covers constrained convex minimization (with $F = \nabla f$), saddle-point problems, Nash equilibria of convex games and complementarity problems, which is why projection methods for it are a standard topic in optimization and numerical analysis.

The classical methods for a monotone, $L$-Lipschitz field (Korpelevich's extragradient method, Tseng's forward–backward–forward method, the subgradient extragradient method of Censor, Gibali and Reich) use a constant step $\lambda < 1/L$ or a backtracking search for one. Yu. Malitsky's paper (arXiv:1502.04968v1; SIAM J. Optim. 25(1), 2015, doi:10.1137/14097238X) introduces the projected reflected gradient method, which needs one projection and one evaluation of $F$ per iteration. Its constant-step version requires $\lambda < (\sqrt2 - 1)/L$, and so requires knowing $L$, which is rarely available in practice (§4, p. 8). Section 4 removes that requirement: Algorithm 4.2 chooses its step from local information only and corrects it when a computable test fails. This mission formalizes the paper's convergence theorem for that algorithm, Theorem 4.4.

Setting

Let $H$ be a real Hilbert space with inner product $\langle\cdot,\cdot\rangle$, $C \subseteq H$ nonempty, closed and convex, and $F : H \to H$. Write $P_C$ for the metric projection onto $C$. The variational inequality (1.1) asks for $x^* \in C$ with

$$\langle F(x^*), x - x^*\rangle \ge 0 \quad\text{for all } x \in C,$$

and $S$ denotes its solution set. The standing conditions (p. 1) are (C1) $S \ne \emptyset$; (C2) $F$ is monotone, $\langle F(x) - F(y), x - y\rangle \ge 0$; (C3) $F$ is $L$-Lipschitz, $\|F(x) - F(y)\| \le L\|x - y\|$, with $L > 0$.

Algorithm 4.2 (pp. 8–9) has parameters $\alpha \in (0, \sqrt2 - 1)$, $\lambda_{-1} > 0$ and a cap $\bar\lambda > 0$, and adopts the convention $0/0 = +\infty$. It starts from $x_0 \in C$ with $y_0 = P_C(x_0 - \lambda_{-1}F(x_0))$, $\lambda_0 = \min\{\alpha\|x_0 - y_0\|/\|F(x_0) - F(y_0)\|, \bar\lambda\}$, $\tau_0 = 1$ and $x_1 = P_C(x_0 - \lambda_0F(y_0))$. For $n \ge 1$ it uses the step rule

$$\lambda(y,\tau) = \min\Big\{\frac{\alpha\|y - y_{n-1}\|}{\|F(y) - F(y_{n-1})\|},\ \frac{1+\tau_{n-1}}{\tau}\lambda_{n-1},\ \bar\lambda\Big\}, \qquad (4.1)$$

forms the trial values $y_n = 2x_n - x_{n-1}$, $\lambda_n = \lambda(y_n, 1)$, $x_{n+1} = P_C(x_n - \lambda_nF(y_n))$, and evaluates a test quantity $t_n$ built from these. If $t_n \le 0$ the trial is accepted. Otherwise, in step (4.i) ($\lambda_n \ge \lambda_{n-1}$) the step is replaced by some $\lambda'_n \in [\lambda_{n-1}, \lambda_n]$ satisfying (4.2), and in step (4.ii) ($\lambda_n < \lambda_{n-1}$) the reflection is replaced by $y'_n = x_n + \tau'_n(x_n - x_{n-1})$ for some $\tau'_n \in (0,1]$ satisfying (4.3), with a step $\lambda'_n$ satisfying (4.4); $x_{n+1}$ is then recomputed. The values $y_n, \lambda_n, \tau_n$ kept after iteration $n$ are the corrected ones. A sequence produced this way is a run of Algorithm 4.2. The Lipschitz constant $L$ appears nowhere in the algorithm.

Formalization targets

Goal: Theorem 4.4 (p. 12)

Under (C1)–(C3), for every run of Algorithm 4.2 there is $x^* \in S$ with

$$x_n \rightharpoonup x^* \qquad\text{and}\qquad y_n \rightharpoonup x^*,$$

where $\rightharpoonup$ is weak convergence in $H$.

Milestones

In the order the proof uses them:

1. Lemma 4.1 (p. 9): step (4.i) is well-defined.
2. Lemma 4.2 (p. 9): step (4.ii) is well-defined.
3. The bound $\lambda(y,\tau) \ge \alpha/L$ from the proof of Lemma 4.2 (p. 9), given $\lambda_{n-1} \ge \alpha/L$ and $\bar\lambda \ge \alpha/L$.
4. Along every run, $0 < \lambda_n \le \bar\lambda$ and $\tau_n \in (0,1]$ (p. 12).
5. Lemma 2.4 (p. 3): $P_M(\bar x + \lambda(x - \bar x)) = \bar x$ for $\bar x = P_Mx$ and $\lambda > 0$.
6. Lemma 4.3, the energy inequality (4.5) for $z \in S$ (p. 10).
7. The bound $\tau_n\lambda_n \le (1 + \tau_{n-1})\lambda_{n-1}$ and the decrease (4.10) (p. 12).
8. Boundedness of $(x_n)$ and $\|x_{n+1} - y_n\| \to 0$, $\|x_{n+1} - x_n\| \to 0$ (p. 12).

The supporting lemmas 2.2 (Minty), 2.3 (Opial) and 2.7 are included as plain items, together with a statement that runs exist from every $x_0 \in C$.

Significance

The result shows that the projected reflected gradient method keeps its weak convergence guarantee when the step is chosen adaptively, without knowledge of $L$, at a cost of at most two projections per iteration (one in most iterations, per the experiments of §5). Adaptive step rules of this kind, based on the local ratio $\|y - y'\|/\|F(y) - F(y')\|$, recur in later work on reflected and golden-ratio methods for variational inequalities (e.g. Malitsky 2020, arXiv:1803.08832), so a formal treatment of the run predicate and its energy inequality is reusable beyond this paper.

The theorem is proved in the paper; it has no machine-checked proof that this mission knows of. The work here is to formalize the definitions of Algorithm 4.2, its well-definedness, the energy inequality and the closing weak-limit argument in a general Hilbert space. The proof as printed has a gap, described under Difficulty; a formal proof would either close it or show that the statement needs an extra hypothesis.

Difficulty

Lemma 4.3 is the bookkeeping core. It has five cases depending on how $\lambda_n$ compares with $\lambda_{n-1}$ and on the sign of a norm expression. In the case where the reflection is changed, the previous iterate must be reinterpreted as the projection of a point on a segment (via Lemma 2.4) before the estimate of the constant-step method applies.

The closing argument needs the step sizes to stay away from zero, so that the limit can be taken in the optimality inequality for the iterates. The paper asserts that $\alpha/L$ is a lower bound of $(\lambda_n)$, "from Lemma 4.2" (p. 12). The proof of Lemma 4.2 shows only that the step rule $\lambda(y,\tau)$ is at least $\alpha/L$, and that only under $\bar\lambda \ge \alpha/L$. In step (4.ii) the algorithm allows any $\tau'_n \in (0,1]$ satisfying (4.3) and any $\lambda'_n \ge \tau'_n\lambda_{n-1}$, so an admissible run can take $\lambda_n$ arbitrarily small at such steps. The step regrows by the factor $1 + \tau_{n-1}$ at later steps, but the printed proof does not control this. The obvious first attempt, citing a uniform lower bound $\lambda_n \ge \alpha/L$, does not work as stated.

Formalization scope

All objects live in one definitions item. $H$ is a general real Hilbert space, never $\mathbb R^n$. The projection $P_C$ is a nearest point chosen by choice, with a junk value when none exists. Every theorem that uses it assumes $C$ closed and convex together with a nonemptiness source, so the junk branch is never reached. Weak convergence is $\langle x_n, v\rangle \to \langle x^*, v\rangle$ for all $v$. (C3) uses a real constant $L > 0$. The convention $0/0 = +\infty$ is encoded by omitting the ratio term from the minimum when its denominator is zero, never by Lean's $x/0 = 0$.

A run is a relation: the corrections $\lambda'_n$, $\tau'_n$ range over all admissible values, and nothing selects the largest. The stopping rule of step 3 is not encoded. A run whose residual vanishes continues as a constant sequence through the branch $t_n \le 0$. Indices $n - 1$ are used only for $n \ge 1$, and $\lambda_{-1}$ is a separate parameter.

The goal assumes exactly (C1)–(C3), $\alpha \in (0, \sqrt2 - 1)$, $\lambda_{-1} > 0$ and $\bar\lambda > 0$. It assumes neither $\bar\lambda \ge \alpha/L$ nor any lower bound on $(\lambda_n)$. Adding such a bound would assume the unproved step of the paper's argument. The milestone on $\lambda(y,\tau) \ge \alpha/L$ carries $\lambda_{n-1} \ge \alpha/L$ and $\bar\lambda \ge \alpha/L$ as explicit, disclosed hypotheses. Lemma 4.3 and (4.10) are stated without (C1) and (C3), which their proofs do not use. A formalization that concluded norm convergence, specialized $H$ to finite dimension, or fixed one particular choice of corrections would be a different theorem.

The development needs the existence and characterization of the metric projection onto a closed convex set, Minty's lemma, Opial's inequality and weak sequential compactness of bounded sets in Hilbert space. These pieces are reusable for other projection methods. Contributions toward any milestone are welcome, as is a resolution of the lower-bound question for $(\lambda_n)$.

Selected references

• Yu. Malitsky, Projected Reflected Gradient Methods for Monotone Variational Inequalities, SIAM J. Optim. 25(1) (2015); read as arXiv:1502.04968v1. https://arxiv.org/abs/1502.04968v1, https://doi.org/10.1137/14097238X
• G. M. Korpelevich, The extragradient method for finding saddle points and other problems, Ekonomika i Matematicheskie Metody 12 (1976), 747–756.
• P. Tseng, A modified forward-backward splitting method for maximal monotone mappings, SIAM J. Control Optim. 38(2) (2000), 431–446. https://doi.org/10.1137/S0363012998338806
• Y. Censor, A. Gibali, S. Reich, The subgradient extragradient method for solving variational inequalities in Hilbert space, J. Optim. Theory Appl. 148 (2011), 318–335. https://doi.org/10.1007/s10957-010-9757-3
• Yu. Malitsky, Golden ratio algorithms for variational inequalities, Math. Program. 184 (2020), 383–410. https://arxiv.org/abs/1803.08832
• Z. Opial, Weak convergence of the sequence of successive approximations for nonexpansive mappings, Bull. Amer. Math. Soc. 73 (1967), 591–597. https://doi.org/10.1090/S0002-9904-1967-11761-0