Convex Optimization

Motivation

Many large convex problems in signal processing, networking and statistics have a separable objective coupled by linear constraints: $K$ blocks of variables, each with its own nonsmooth regularizer and its own feasible set, tied together by $\sum_k E_kx_k=q$. The alternating direction method of multipliers (ADMM) is the standard tool for such problems when $K=2$. For $K\ge3$ the direct extension of ADMM can diverge, as shown by Chen, He, Ye and Yuan (2016), and convergence guarantees for multi-block variants typically require strong convexity or extra correction steps.

Hong, Chang, Wang, Razaviyayn, Ma and Luo propose the block successive upper-bound minimization method of multipliers (BSUM-M), which updates each primal block by minimizing a local upper bound of the augmented Lagrangian, followed by a dual gradient step, and prove that it converges for any number of blocks without strong convexity of the objective. They also analyse a randomized version, RBSUM-M, in which each iteration updates a single randomly chosen block, primal or dual. Randomized block selection matters when data or computation are distributed and not every block is available at every step. This mission formalizes the randomized half of their main theorem; the cyclic half is a separate mission.

Setting

Let $x=(x_1,\dots,x_K)$ with $x_k\in\mathbb R^{n_k}$, and consider

$$\min_x\ f(x)=g(x)+\sum_{k=1}^K h_k(x_k)\quad\text{s.t.}\quad Ex=\sum_{k=1}^K E_kx_k=q,\quad x_k\in X_k. \tag{1.1}$$

The standing assumptions (Assumption A) are: $g(x)=\ell(Ax)+\langle x,b\rangle$ with $\ell$ strictly convex and continuously differentiable; $h_k(x_k)=\lambda_k\|x_k\|_1+\sum_J w_J\|x_{k,J}\|_2$ with nonnegative weights; each $X_k=\{x_k\mid C_kx_k\le c_k\}$ is a compact polyhedron; the problem is feasible and its primal and dual optimal values are attained.

For $\rho>0$ the augmented Lagrangian and augmented dual function are

$$L(x;y)=f(x)+\langle y,q-Ex\rangle+\frac\rho2\|q-Ex\|^2,\qquad d(y)=\min_{x\in X}L(x;y),$$

and $X(y)$ is the set of minimizers of $L(\cdot;y)$ over $X=\prod_kX_k$. The method works with approximation functions $u_k(v_k;x)$ (Assumption B): $u_k(\cdot;x)$ majorizes $g+\frac\rho2\|E\cdot-q\|^2$ in block $k$, touches it with matching gradient at $x_k$, is strongly convex with modulus $\gamma_k$ and has an $L_k$-Lipschitz gradient.

RBSUM-M. Fix probabilities $p_0,\dots,p_K>0$ summing to 1. At iteration $t\ge1$ draw $k$ with probability $p_k$. If $k=0$, take a dual step $y^{t+1}=y^t+\alpha^t(q-Ex^t)$; otherwise replace block $k$ by

$$x_k^{t+1}=\arg\min_{x_k\in X_k}\ u_k(x_k;x^t)-\langle y^t,E_kx_k\rangle+h_k(x_k),$$

leaving the other blocks and $y$ unchanged. The analysis uses the block steps $\hat x_k^{t+1}$ (the minimizer above, computed for every $k$), the dual step $\hat y^{t+1}=y^t+\alpha^t(q-Ex^t)$, the gaps $\Delta_d^t=d^*-d(y^t)$, $\Delta_p^t=L(x^t;y^t)-d(y^t)$, and the proximal gradient $\tilde\nabla_xL(x;y)=x-\operatorname{prox}_{h+\iota_X}\big(x-\nabla_x(L(x;y)-h(x))\big)$. The error bound with constant $\tau$ is $\operatorname{dist}(x,X(y))\le\tau\|\tilde\nabla_xL(x;y)\|$ for all $y$ and $x\in X$.

Formalization targets

Goal: Theorem 2.1, part 2

Under Assumptions A and B and the error bound, if the stepsizes are either constant and sufficiently small, or satisfy $\sum_t\alpha^t=\infty$, $\alpha^t\to0$, then with probability 1

$$\|Ex^t-q\|\to0,\qquad\|x^t-x^{t+1}\|\to0,\qquad\|x^t-\bar x^t\|\to0,$$

where $\bar x^t$ is the point of $X(\hat y^t)$ nearest to $x^t$, and every limit point of $(x^t,y^t)$ is a primal and dual optimal pair. The threshold for "sufficiently small" is existential and depends only on the problem data, the $u_k$, the error-bound constant and $p$.

Milestones

In the order of the paper: Lemma 2.1 (differentiability of $d$, $\nabla d(y)=q-Ex(y)$, the $1/\rho$-Lipschitz bound); Lemma 2.3(2) (expected decrease of $L$); Lemma 2.4(2) (proximal gradient bounded by the block and dual steps); Lemmas 2.5(2) and 2.6(2) (expected change of the dual and primal gaps); and from the proof of Theorem 2.1, (2.24)–(2.25) (expected change of $\Delta_p+\Delta_d$), (2.28) (constraint violation), (2.29) (the descent estimate with explicit coefficients) and (2.40) (one-step change of $\|E\bar x^t-q\|$).

Significance

The result shows that a randomized, single-block-per-iteration primal-dual method converges almost surely to primal-dual optimal solutions of a multi-block linearly constrained problem, under assumptions that allow a rank-deficient $A$, nonsmooth group-sparse regularizers, and inexact (majorized) block subproblems. The potential function $\Delta_p+\Delta_d$ and the way the error bound converts it into a supermartingale estimate are reused in later analyses of multi-block and randomized ADMM-type methods.

The result is proved in the paper; to our knowledge it has no machine-checked proof. Formalizing it produces a checked convergence proof for a randomized augmented-Lagrangian method, a reusable formal treatment of the augmented dual function of a convex program with compact polyhedral constraints, and a precise record of the conventions under which the published statements hold (see the scope section).

Difficulty

The obvious route is to show that $L(x^t;y^t)$ decreases. It does not: a dual step increases $L$ by $\alpha^t\|q-Ex^t\|^2$, and without strong convexity there is no direct control of the distance to the solution set. The argument has to combine the primal and dual gaps into one potential, bound the residual term $\|Ex^t-E\bar x^{t+1}\|$ through an error bound that holds without strong convexity, and then pass from a conditional-expectation inequality to almost-sure convergence. Under diminishing stepsizes the supermartingale argument only yields $\liminf\|E\bar x^t-q\|=0$, and upgrading this to a limit needs a separate pathwise argument on how fast $\nabla d(\hat y^t)$ can move.

Formalization scope

Blocks are EuclideanSpace ℝ (Fin (n k)) and the full vector lives in PiLp 2, so norms are Euclidean. $\ell$ is real-valued on all of $\mathbb R^p$. The augmented dual function is $d(y)=\min_{x\in X}L(x;y)$ (the page's (1.9) prints an unconstrained minimum of $g$ plus penalty; the constrained reading is the one Lemma 2.1 and (2.19) use). The prox in the proximal gradient includes the constraint set $X$, which is what the proof of Lemma 2.4 needs; the error bound of Lemma 2.2 is a hypothesis of the goal, in its global form on $X$, and is not derived. Lemma 2.2 itself is not part of the mission.

Conditional expectations given $z^t=(x^t,y^t)$ are written as the explicit average $p_0F(x,\hat y)+\sum_kp_kF((\hat x_k,x_{-k}),y)$ over the random index, at every state; this is $\mathbb E[F(z^{t+1})\mid z^t]$ because the index drawn at iteration $t$ is independent of $z^t$. The goal is stated on an arbitrary probability space carrying mutually independent indices with law $p$, for every run from a deterministic start $x^1\in X$, $y^1$. Iterations are indexed from $t=1$ as on the page.

Two printed statements are corrected and labelled in their items: Lemma 2.1 asserts $Ax$ (not each $A_kx_k$) constant over $X(y)$, and Lemma 2.6(2) uses $\alpha^t$ where the page prints $\alpha^r$. Statements are not trivialized: the setting is satisfiable (for instance $K=1$, $g\equiv0$, $h=0$, $E=0$, $X=[-1,1]$, $u(v;x)=(v-x_1)^2/2$), the stepsize threshold precedes the probability space and the run, and the solution sets entering distances are nonempty under Assumption A. A formalization in which the error bound or the gap bounds are assumed in a form that already contains the conclusion would not count.

A complete development needs: Danskin-type differentiability of $d$ and the $1/\rho$-smoothness of the augmented dual; first-order optimality and nonexpansiveness of the prox of $h+\iota_X$; and the Robbins–Siegmund almost-supermartingale convergence theorem in Mathlib's conditional expectation framework, together with a bridge from the explicit index average to condExp. The last two are reusable well beyond this mission. Contributions to any of the milestones, or to these general lemmas, are welcome.

Selected references

• M. Hong, T.-H. Chang, X. Wang, M. Razaviyayn, S. Ma, Z.-Q. Luo, A Block Successive Upper Bound Minimization Method of Multipliers for Linearly Constrained Convex Optimization, arXiv:1401.7079v1, 2014; Mathematics of Operations Research 45(3), 2020. https://arxiv.org/abs/1401.7079, https://doi.org/10.1287/moor.2019.1010
• M. Hong, Z.-Q. Luo, On the linear convergence of the alternating direction method of multipliers, Mathematical Programming 162, 2017. https://doi.org/10.1007/s10107-016-1034-2
• C. Chen, B. He, Y. Ye, X. Yuan, The direct extension of ADMM for multi-block convex minimization problems is not necessarily convergent, Mathematical Programming 155, 2016. https://doi.org/10.1007/s10107-014-0826-5
• H. Robbins, D. Siegmund, A convergence theorem for non negative almost supermartingales and some applications, in Optimizing Methods in Statistics, Academic Press, 1971. https://doi.org/10.1016/B978-0-12-604550-5.50015-8
• M. Razaviyayn, M. Hong, Z.-Q. Luo, A unified convergence analysis of block successive minimization methods for nonsmooth optimization, SIAM Journal on Optimization 23(2), 2013. https://doi.org/10.1137/120891009

Theorem numbers, equation numbers and page numbers in this mission refer to arXiv:1401.7079v1 (printed page = PDF page − 1); the journal version renumbers them.

Motivation

Many problems in signal processing, machine learning and power systems have the form of a convex objective that is a sum of a smooth coupled term and nonsmooth block-separable terms, minimized subject to linear equality constraints that couple $K$ blocks of variables: basis pursuit, demand response control in smart grids, and distributed estimation are the examples of Hong, Chang, Wang, Razaviyayn, Ma, Luo. The alternating direction method of multipliers (ADMM) handles two blocks well, but for $K\ge3$ blocks the directly extended Gauss–Seidel ADMM can diverge (Chen, He, Ye, Yuan, 2016), and convergence proofs for multi-block schemes usually need strong convexity or extra correction steps.

The block successive upper-bound minimization method of multipliers (BSUM-M) is the paper's answer: each block is updated by minimizing a local upper bound of the augmented Lagrangian, which makes the subproblems simple (for instance closed-form proximal steps), and the multiplier is updated by a dual gradient step with a small or diminishing stepsize. The main theorem shows that this scheme converges to primal and dual optimal solutions for any number of blocks without strong convexity, under an error bound.

Setting

The variable is $x=(x_1,\dots,x_K)$, $x_k\in\mathbb R^{n_k}$, with the Euclidean norm $\|x\|^2=\sum_k\|x_k\|^2$. Problem (1.1) is

$$\min_x\ f(x)=g(x)+\sum_{k=1}^K h_k(x_k)\quad\text{s.t.}\quad Ex=\sum_{k=1}^KE_kx_k=q,\quad x_k\in X_k .$$

Under the standing Assumption A, $g(x)=\ell(Ax)+\langle x,b\rangle$ with $\ell$ strictly convex and continuously differentiable and $A$ any matrix; $h_k(x_k)=\lambda_k\|x_k\|_1+\sum_Jw_J\|x_{k,J}\|_2$ is a mixed $\ell_1/\ell_2$ norm with nonnegative weights; each $X_k=\{x_k\mid C_kx_k\le c_k\}$ is a compact polyhedron and $X=\prod_kX_k$; the problem is feasible and the dual optimal value is attained.

For $\rho>0$, the augmented Lagrangian is $L(x;y)=f(x)+\langle y,q-Ex\rangle+\frac\rho2\|q-Ex\|^2$, the augmented dual is $d(y)=\min_{x\in X}L(x;y)$, and $X(y)=\arg\min_{x\in X}L(x;y)$. For $x\in X$, $\bar x$ denotes the point of $X(y)$ nearest to $x$.

Assumption B concerns the approximation functions $u_k(v_k;x)$: each upper-bounds $G(x)=g(x)+\frac\rho2\|Ex-q\|^2$ along block $k$, is tight with matching gradient at $v_k=x_k$, is continuous, strongly convex in $v_k$ with a uniform modulus $\gamma_k>0$, and has an $L_k$-Lipschitz gradient in $v_k$.

BSUM-M (1.12). From $x^1\in X$ and any $y^1$, for $r\ge1$:

$$y^{r+1}=y^r+\alpha^r(q-Ex^r),\qquad x_k^{r+1}=\arg\min_{x_k\in X_k}u_k(x_k;w_k^{r+1})-\langle y^{r+1},E_kx_k\rangle+h_k(x_k),$$

for $k=1,\dots,K$ in order, with the Gauss–Seidel point $w_k^{r+1}=(x_1^{r+1},\dots,x_{k-1}^{r+1},x_k^r,\dots,x_K^r)$.

The proximal gradient is $\tilde\nabla_xL(x;y)=x-\operatorname{prox}(x-\nabla_x(L(x;y)-h(x)))$ with $\operatorname{prox}(z)=\arg\min_{u\in X}h(u)+\frac12\|z-u\|^2$. The error bound asks for $\tau>0$ with $\operatorname{dist}(x,X(y))\le\tau\|\tilde\nabla_xL(x;y)\|$ for all $x\in X$ and all $y$.

Formalization targets

Goal: Theorem 2.1, part 1

Under Assumptions A and B and the error bound, there is $\bar\alpha>0$ such that for stepsizes $\alpha^r>0$ that are either constant with $\alpha^r=\alpha\le\bar\alpha$, or satisfy $\sum_r\alpha^r=\infty$ and $\alpha^r\to0$, every BSUM-M run satisfies

$$\lim_{r\to\infty}\|Ex^r-q\|=0,\qquad\lim_{r\to\infty}\|x^r-x^{r+1}\|=0,\qquad\lim_{r\to\infty}\|x^r-\bar x^r\|=0,$$

and every limit point of $\{(x^r,y^r)\}$ is a primal and dual optimal solution. The threshold $\bar\alpha$ is not specified; no rate is claimed.

Milestones

Lemma 2.3(1) (one sweep decreases $L(\cdot;y^{r+1})$ by $\gamma\|x^r-x^{r+1}\|^2$), Lemma 2.5(1) (dual gap decrease), Lemma 2.6(1) (primal gap bound), Lemma 2.4(1) (the proximal gradient at $x^r$ is $O(\|x^{r+1}-x^r\|)$) and Lemma 2.1 (differentiability of $d$, $\nabla d(y)=q-Ex(y)$, and Lipschitz continuity of $\nabla d$ on superlevel sets).

Significance

The theorem gives a convergent multi-block method of multipliers whose block subproblems may be replaced by any majorizer satisfying Assumption B, such as a linearized proximal step. It covers the cyclic Gauss–Seidel order, where the direct multi-block ADMM can fail, and it does not require strong convexity of $f$ (the matrix $A$ need not have full column rank). The paper's randomized variant (part 2 of the same theorem) is a separate mission.

The result is proved in the paper (the proof of part 1 is described as following the steps written out for part 2). It has not been machine-checked. Formalization adds a verified account of the potential-function argument for a primal-dual method with inexact block updates, and checks three printed statements that need correction (see the scope section): a formal proof settles them.

Difficulty

The obvious approach is to treat BSUM-M as an inexact dual gradient ascent on $d$, but a single Gauss–Seidel sweep does not compute $x(y^{r+1})\in X(y^{r+1})$, so the dual step uses a gradient at the wrong point. Controlling that error requires relating the primal step length to the distance from $X(y)$, which fails without an error bound: $L(\cdot;y)$ is not strongly convex, $X(y)$ need not be a singleton, and the proximal gradient can be small far from $X(y)$. The cyclic order adds a second gap: block $k$ is linearized at $w_k^{r+1}$, not at $x^r$. With a constant stepsize the dual ascent must not outrun the primal descent, which is where the smallness threshold enters.

Formalization scope

Citation basis: the arXiv preprint arXiv:1401.7079v1 (28 Jan 2014); every index refers to that version, not to the revised journal version (Mathematics of Operations Research, 2020).

Representation: blocks are indexed by Fin K with $K\ge1$; the block space is PiLp 2 of Euclidean spaces (the Euclidean norm, not the sup norm). $\ell$ is real valued and $C^1$ on all of $\mathbb R^p$, which specializes the paper's extended-valued $\ell$. $d$, $d^*$ and $f^*$ are real infima and suprema over nonempty compact sets, where they are attained. The algorithm is a predicate on sequences indexed from $r=1$; argmin steps are minimality conditions, the minimizer being unique. The initial point is required to lie in $X$. "Sufficiently small" is an existential threshold quantified before the stepsizes and the run. $\|x^r-\bar x^r\|$ is the distance from $x^r$ to $X(y^r)$. Limit points are cluster points of the joint sequence; boundedness of $\{y^r\}$ is neither assumed nor concluded.

Interpretations and corrections, each labelled in the item's statement:

• the proximity operator in (2.3) includes the constraint set $X$, as the paper's optimality relation (2.14) requires;
• $d(y)$ is $\min_{x\in X}L(x;y)$; (1.9) prints $g(x)$ and an unconstrained minimum;
• (2.17) has $Ex^r$ where its proof (2.19) has $Ex^{r-1}$; the proved inequality is stated;
• (2.12) evaluates the proximal gradient at $y^r$ and is false at $r=1$; it is stated at $y^{r+1}$;
• Lemma 2.1's clause "$A_kx_k$ constant on $X(y)$" fails for non-separable $\ell$; it is stated as "$Ax$ constant".

The error bound is a hypothesis of the goal, as in the paper, and is not derived from Assumption A: the paper's Lemma 2.2 does not hold under Assumption A alone. A formalization that makes the hypotheses unsatisfiable, drops the run's dependence on $y^{r+1}$, or lets $\bar\alpha$ depend on the run is not this theorem. A sanity check that the hypotheses are satisfiable (a one-block instance with a constant run) is part of the drafting record.

Infrastructure needed: proximity operators of convex functions restricted to polyhedra and their nonexpansiveness; Danskin-type differentiability of a parametric minimum over a compact set; smoothness of the augmented dual; elementary convergence lemmas for nonnegative sequences with summable decrements. These pieces are reusable for other augmented Lagrangian and ADMM analyses. Proofs of the milestones are welcome independently of the goal.

Selected references

• M. Hong, T.-H. Chang, X. Wang, M. Razaviyayn, S. Ma, Z.-Q. Luo, A Block Successive Upper Bound Minimization Method of Multipliers for Linearly Constrained Convex Optimization, arXiv:1401.7079v1, 2014; Mathematics of Operations Research 45(3), 2020. https://arxiv.org/abs/1401.7079 , https://doi.org/10.1287/moor.2019.1010
• M. Hong, Z.-Q. Luo, On the linear convergence of the alternating direction method of multipliers, Mathematical Programming 162, 2017. https://doi.org/10.1007/s10107-016-1034-2
• M. Razaviyayn, M. Hong, Z.-Q. Luo, A unified convergence analysis of block successive minimization methods for nonsmooth optimization, SIAM Journal on Optimization 23(2), 2013. https://doi.org/10.1137/120891009
• C. Chen, B. He, Y. Ye, X. Yuan, The direct extension of ADMM for multi-block convex minimization problems is not necessarily convergent, Mathematical Programming 155, 2016. https://doi.org/10.1007/s10107-014-0826-5

Motivation

Best-subset selection with ridge shrinkage, $\min_\beta \tfrac12\|y-X\beta\|_2^2+\lambda_0\|\beta\|_0+\lambda_2\|\beta\|_2^2$, is a mixed-integer program that statisticians want to solve to certified optimality at the scale of modern data ($p$ up to $10^7$ features). Hazimeh, Mazumder and Saab (arXiv:2004.06152, Mathematical Programming 2022) built a branch-and-bound solver, L0BnB, whose node relaxations are solved in the primal space by active-set coordinate descent instead of by an interior-point method. Branch-and-bound prunes a node only with a valid dual bound, a certified lower bound on the node's relaxation value. A primal method returns an approximate minimizer, not such a bound, so the solver must turn an inexact primal point into a dual feasible point and needs to know how much is lost in doing so. This mission formalizes the paper's answer, its main theorem (Theorem 3).

Setting

Data are $X\in\mathbb R^{n\times p}$ with columns $X_1,\dots,X_p$, $y\in\mathbb R^n$, and parameters $\lambda_0,\lambda_2,M>0$. The reverse Huber penalty is $\mathcal B(t)=|t|$ for $|t|\le1$ and $(t^2+1)/2$ for $|t|\ge1$. The penalty $\psi(b)$ equals $\psi_1(b)=2\lambda_0\mathcal B(b\sqrt{\lambda_2/\lambda_0})$ if $\sqrt{\lambda_0/\lambda_2}\le M$ and $\psi_2(b)=(\lambda_0/M+\lambda_2M)|b|$ otherwise. The reduced relaxation (5) is

$$\min_{\beta\in\mathbb R^p}\ F(\beta)=\tfrac12\|y-X\beta\|_2^2+\sum_{i}\psi(\beta_i)\quad\text{s.t.}\quad \|\beta\|_\infty\le M .$$

Throughout, the columns of $X$ and $y$ have unit $\ell_2$ norm (the standing assumption of Section 3 of the paper).

Algorithm 2 (active-set coordinate descent) returns a point $\hat\beta$ in the box such that the set

$$V=\{i\notin\operatorname{Supp}(\hat\beta):\ 0\ne\operatorname*{arg\,min}_{|t|\le M}F(\hat\beta_1,\dots,t,\dots,\hat\beta_p)\}$$

is empty: no coordinate outside the support wants to move. Write $\hat r=y-X\hat\beta$, $k=\|\hat\beta\|_0$, and let $\beta^*$ be an optimal solution of (5) with $r^*=y-X\beta^*$. The primal gap is $\epsilon=\|X(\beta^*-\hat\beta)\|_2$.

The two duals of (5) (Theorem 2 of the paper) are, with $v(\alpha,\gamma_i)=[(\alpha^\top X_i-\gamma_i)^2/(4\lambda_2)-\lambda_0]_++M|\gamma_i|$,

$$h_1(\alpha,\gamma)=-\tfrac12\|\alpha\|_2^2-\alpha^\top y-\textstyle\sum_i v(\alpha,\gamma_i),\qquad h_2(\rho,\mu)=-\tfrac12\|\rho\|_2^2-\rho^\top y-M\|\mu\|_1,$$

the latter under $|\rho^\top X_i|-\mu_i\le\lambda_0/M+\lambda_2M$. The dual variables built from $\beta^*$ are $\alpha^*=\rho^*=-r^*$, $\gamma^*_i=\mathbb 1_{[|\beta^*_i|=M]}(\alpha^{*\top}X_i-2M\lambda_2\operatorname{sign}(\alpha^{*\top}X_i))$ and $\mu^*_i=\mathbb 1_{[|\beta^*_i|=M]}(|\rho^{*\top}X_i|-\lambda_0/M-\lambda_2M)$. The dual points built from $\hat\beta$ are $\hat\alpha=\hat\rho=-\hat r$, with $\hat\gamma$ a maximizer of $h_1(\hat\alpha,\cdot)$ (25) and $\hat\mu$ a maximizer of $h_2(\hat\rho,\cdot)$ under the constraints (27).

Formalization targets

Goal: Theorem 3 with the proof's constants

If $\sqrt{\lambda_0/\lambda_2}\le M$, with $c_i=(2\lambda_2)^{-1}$ when $|\beta^*_i|<M$ and $c_i=M$ when $|\beta^*_i|=M$,

$$h_1(\hat\alpha,\hat\gamma)\ \ge\ h_1(\alpha^*,\gamma^*)-2\epsilon-\tfrac12\epsilon^2-\sum_{i\in\operatorname{Supp}(\hat\beta)}\big(c_i\epsilon+(4\lambda_2)^{-1}\epsilon^2\big).\qquad(55)$$

If $\sqrt{\lambda_0/\lambda_2}>M$,

$$h_2(\hat\rho,\hat\mu)\ \ge\ h_2(\rho^*,\mu^*)-\epsilon(2+Mk)-\tfrac12\epsilon^2.\qquad(59)$$

The paper states these as $-kO(\epsilon)-kO(\epsilon^2)$ (29) and $-kO(\epsilon)-O(\epsilon^2)$ (30); the goal states the expressions its proof establishes.

Milestones

In the order the proof uses them: the closed-form coordinate updates (14) and (15); Proposition 3, $V=\{i\notin\operatorname{Supp}(\hat\beta): |\langle\hat r,X_i\rangle|>c(\lambda_0,\lambda_2,M)\}$; the bound $\|\alpha^*\|_2\le1$; Lemma 2, which bounds $v(\hat\alpha,\hat\gamma_i)$ coordinate by coordinate and makes it vanish off the support; inequality (53); the closed form (28) of $\hat\mu$; the optimality conditions (58); the bound (57) $|\hat\mu_i|\le\epsilon+|\mu^*_i|$ on the support; and inequality (56).

Significance

The result. In both regimes the loss of the dual bound is controlled by the primal gap times the sparsity $k$ of the iterate, with constants depending only on $M$ and $\lambda_2$; the number of features $p$ does not enter. Since L0BnB seeks solutions with $k\ll p$, the cheap dual bound obtained from an inexact primal solution is nearly as good as the exact one, which is what allows pruning without an interior-point solve at every node. The paper notes that with plain coordinate descent in place of Algorithm 2 the same argument gives $p$ in place of $k$.

Formalizing it. The result is proved in the paper; it has no machine-checked proof. The mission produces a checked version of the main theorem with explicit constants, a formal model of what "output of an active-set method" means for the analysis, and reusable closed forms for the boxed soft-thresholding updates of an $\ell_1$- or reverse-Huber-penalized box-constrained least squares problem.

Difficulty

The bound is not a consequence of weak duality alone: weak duality says only that each dual value is below the primal optimum, and says nothing about how far the constructed dual point is from the dual optimum. The obvious estimate, a Lipschitz bound on $h_1$ or $h_2$ summed over all coordinates, gives a loss proportional to $p$. Getting $k$ instead requires showing that the dual contribution of every coordinate outside $\operatorname{Supp}(\hat\beta)$ vanishes exactly, which uses the emptiness of $V$ through Proposition 3 and the closed forms (14)–(15) of one-dimensional nonsmooth box-constrained problems. On the support, the loss must be bounded using the structure of $\gamma^*$ and $\mu^*$ at coordinates where $\beta^*$ hits the box, which is where $c_i=M$ and (58) enter.

Formalization scope

Indices are Fin p; norms are explicit sums; $\operatorname{sign}$ is Real.sign with $\operatorname{sign}(0)=0$; $[a]_+$ is max a 0. The unit norms of the columns of $X$ and of $y$ are hypotheses, never built into the definitions, and each theorem assumes only the ones it uses. $\beta^*$ is any point of the box minimizing $F$ over the box. The output of Algorithm 2 is modelled by its two properties used in the paper: $\|\hat\beta\|_\infty\le M$ and $V=\emptyset$ ("$0\ne\arg\min$" encoded as "$0$ is not a minimizer"); the iterations are not formalized. $\hat\gamma$ and $\hat\mu$ are arbitrary maximizers satisfying (25) and (27). The dual variables (23)–(24) are defined by their formulas; their optimality (Theorem 2, a separate mission) is neither assumed nor needed, and the goal is an inequality between explicit numbers.

Instantiations of $O(\cdot)$. (29) is stated as (55): $kO(\epsilon)+kO(\epsilon^2)$ becomes $2\epsilon+\tfrac12\epsilon^2+\sum_{i\in\operatorname{Supp}(\hat\beta)}(c_i\epsilon+(4\lambda_2)^{-1}\epsilon^2)$ with $c_i\in\{(2\lambda_2)^{-1},M\}$ decided by $\beta^*$. (30) is stated as (59): $kO(\epsilon)+O(\epsilon^2)$ becomes $\epsilon(2+Mk)+\tfrac12\epsilon^2$. The proof's final rearrangement of (29), which counts $|\hat\beta_i|$ instead of $|\beta^*_i|$, is not used.

A formalization that takes $\hat\beta=\beta^*$, $\hat\gamma=\gamma^*$ or $\hat\mu=\mu^*$, that drops the hypothesis $V=\emptyset$, or that replaces the constants by an existential "$\exists C$" would be trivial or a different statement; the goal quantifies over every optimal $\beta^*$, every box-feasible $\hat\beta$ with $V=\emptyset$ and every maximizer $\hat\gamma$, $\hat\mu$, with the explicit constants above.

Needed infrastructure: one-dimensional convex minimization on an interval with piecewise penalties, Cauchy–Schwarz for finite sums, and first-order optimality conditions for box-constrained composite problems. Proofs of any milestone are welcome independently.

Selected references

• H. Hazimeh, R. Mazumder, A. Saab, Sparse Regression at Scale: Branch-and-Bound rooted in First-Order Optimization, arXiv:2004.06152v2, 2021; Mathematical Programming 196 (2022). https://arxiv.org/abs/2004.06152
• P. Tseng, Convergence of a block coordinate descent method for nondifferentiable minimization, J. Optim. Theory Appl. 109 (2001). https://doi.org/10.1023/A:1017501703105
• J. Friedman, T. Hastie, R. Tibshirani, Regularization paths for generalized linear models via coordinate descent, J. Stat. Softw. 33 (2010). https://doi.org/10.18637/jss.v033.i01

Motivation

Second order dynamical systems with vanishing damping are continuous-time models of accelerated first-order methods. The system $\ddot x + \frac{\alpha}{t}\dot x + \nabla g(x) = 0$ is the continuous limit of Nesterov's accelerated gradient method (Su, Boyd, Candès 2016), and its trajectories converge weakly to a minimizer of a convex $g$ when $\alpha>3$ (Attouch, Chbani, Peypouquet, Redont 2018). Two modifications of this system have been studied separately. A Hessian-driven damping term $\beta\nabla^2 g(x)\dot x$ damps oscillations and keeps the fast rates (Attouch, Peypouquet, Redont 2016). A Tikhonov regularization term $\epsilon(t)x$ with $\epsilon(t)\to0$ selects one minimizer, the one of minimum norm, and can turn weak convergence into strong convergence (Attouch, Chbani, Riahi 2018).

Boţ, Csetnek and László (arXiv:1911.12845v2, Math. Program. 2021) study the system with both terms. Their results split into two regimes according to how fast $\epsilon$ decays. This mission covers the slow-decay regime of §4.1, where $\int^{+\infty}\epsilon(t)/t\,dt=+\infty$ and the trajectory approaches the minimum-norm minimizer in a weighted average sense.

Timeline. 2016: Su, Boyd and Candès derive the continuous model of Nesterov's method. 2016: Attouch, Peypouquet and Redont add Hessian-driven damping ($\alpha\ge3$, $\beta>0$) and prove fast rates and weak convergence. 2018: Attouch, Chbani and Riahi add Tikhonov regularization without Hessian damping and prove, among other results, strong ergodic convergence to the minimum-norm minimizer when $\int\epsilon(t)/t\,dt=+\infty$. 2020: Boţ, Csetnek and László combine both terms and extend that ergodic result (their Theorem 4.2).

Setting

Let $\mathcal H$ be a real Hilbert space, $t_0>0$, $\alpha>0$, $\beta\ge0$, and $u_0,v_0\in\mathcal H$. The data satisfy the paper's General assumption:

• $g:\mathcal H\to\mathbb R$ is convex and twice Fréchet differentiable, its gradient $\nabla g$ is Lipschitz continuous on bounded sets, and $\operatorname{argmin} g\neq\emptyset$;
• $\epsilon:[t_0,+\infty)\to[0,+\infty)$ is nonincreasing, of class $C^1$, and $\lim_{t\to+\infty}\epsilon(t)=0$.

A global $C^2$-solution of system (5) is a twice continuously differentiable $x:[t_0,+\infty)\to\mathcal H$ with

$$\ddot x(t)+\frac{\alpha}{t}\dot x(t)+\beta\nabla^2 g(x(t))\dot x(t)+\nabla g(x(t))+\epsilon(t)x(t)=0\quad(t\ge t_0),\qquad x(t_0)=u_0,\ \dot x(t_0)=v_0 .$$

The set $\operatorname{argmin} g$ is nonempty, closed and convex, so it has a unique element of minimum norm, the minimum-norm minimizer $x^*=\operatorname{argmin}\{\|x\|:x\in\operatorname{argmin} g\}$. For $\epsilon>0$ the Tikhonov approximation curve is

$$x_\epsilon=\operatorname*{argmin}_{x\in\mathcal H}\Big(g(x)+\frac{\epsilon}{2}\|x\|^2\Big),$$

the unique minimizer of a strongly convex function. Finally, $h_{x^*}(t)=\tfrac12\|x(t)-x^*\|^2$ measures the distance of the trajectory to $x^*$, with derivative $\dot h_{x^*}(t)=\langle\dot x(t),x(t)-x^*\rangle$.

Formalization targets

Goal: Theorem 4.2 (p. 18)

If $\int_{t_0}^{+\infty}\frac{\epsilon(t)}{t}\,dt=+\infty$ and $\alpha>0$, then every global $C^2$-solution satisfies

$$\lim_{t\to+\infty}\frac{1}{\int_{t_0}^{t}\frac{\epsilon(s)}{s}\,ds}\int_{t_0}^{t}\frac{\epsilon(s)}{s}\|x(s)-x^*\|^2\,ds=0 \qquad\text{and}\qquad \liminf_{t\to+\infty}\|x(t)-x^*\|=0 .$$

No rate and no constant is asserted, so the goal does not depend on any particular choice of $\epsilon$.

Milestones

1. Lemma 4.1 (p. 16): for $\alpha>0$, $\beta\ge0$, the velocity is bounded, $\frac1t\|\dot x(t)\|^2\in L^1([t_0,+\infty))$, and $\sup_{t\ge t_0}\frac1t|\dot h_{x^*}(t)|<+\infty$ for every $x^*\in\operatorname{argmin} g$.
2. §4, p. 18: $\|x_\epsilon\|\le\|x^*\|$ for every $\epsilon>0$.
3. §4, p. 18: $\lim_{\epsilon\to0^+}x_\epsilon=x^*$ (stated in the paper as well known).
4. (46) (p. 20): under the hypotheses of Theorem 4.2 there is $C>0$ with

$$\int_{t_0}^{t}\frac{\epsilon(s)}{s}\Big(h_{x^*}(s)-\frac12\big(\|x^*\|^2-\|x_{\epsilon(s)}\|^2\big)\Big)ds\le C\quad\text{for every }t\ge t_0 .$$

Significance

The result shows that slow Tikhonov regularization still selects the minimum-norm minimizer when Hessian damping is added: the weighted time average of $\|x(t)-x^*\|^2$ tends to zero and the trajectory comes arbitrarily close to $x^*$ infinitely often. The hypothesis is only $\alpha>0$, well below the threshold $\alpha\ge3$ that the paper needs for its fast rates, and no growth condition on $\epsilon$ beyond the divergence of $\int\epsilon(t)/t\,dt$ is imposed. It complements Theorem 4.4 of the same paper (a separate mission), which reaches full strong convergence under fast decay, $\int\epsilon(t)/t\,dt<+\infty$, plus further conditions.

On the formal side, the theorem and its milestones are proved in the paper and in the cited literature, but none of them has a machine-checked proof. The work splits into an energy estimate for a nonautonomous second order ODE in a Hilbert space (Lemma 4.1), two facts about the Tikhonov curve that underlie the whole theory of Tikhonov regularization of convex problems (milestones 2 and 3), an integrated differential inequality (46), and a l'Hospital-type averaging step. The Tikhonov-curve facts are reusable for any formal treatment of viscosity selection and minimum-norm solutions.

Difficulty

The obvious attempt, a Lyapunov function that decreases along the trajectory and controls $\|x(t)-x^*\|$, fails because $x^*$ is not a stationary point of the perturbed system: $\nabla g(x^*)+\epsilon(t)x^*=\epsilon(t)x^*\neq0$ in general. The distance to $x^*$ therefore need not decrease, and pointwise convergence is not available under the slow-decay hypothesis alone. The comparison point has to move along the Tikhonov curve $x_{\epsilon(t)}$, whose convergence to $x^*$ comes without a rate. The available bound controls only a weighted integral of $h_{x^*}$ corrected by $\|x^*\|^2-\|x_{\epsilon(s)}\|^2$, so the passage from (46) to the goal needs both $x_{\epsilon(t)}\to x^*$ and the divergence of the weight. In the formal setting the Tikhonov-curve limit is itself a nontrivial weak-compactness argument in a Hilbert space, which the paper does not spell out.

Formalization scope

• $\mathcal H$ is a real Hilbert space (InnerProductSpace ℝ H, CompleteSpace H). $\nabla g$ is gradient g, and $\nabla^2 g(x)v$ is the Fréchet derivative of gradient g at $x$ applied to $v$. "Twice Fréchet differentiable" means $g$ and $\nabla g$ are differentiable; Lipschitz continuity on bounded sets is stated on every closed ball around the origin.
• Trajectories are maps $\mathbb R\to\mathcal H$ with explicit velocity and acceleration maps; derivatives are taken within $[t_0,+\infty)$ (one-sided at $t_0$), the acceleration is continuous there, and values before $t_0$ are irrelevant. The same holds for $\epsilon$ and its derivative.
• Every theorem is stated for every global $C^2$-solution of (5). Existence and uniqueness of the solution is Theorem 2.1 of the paper, a milestone of the first mission of this series, so the hypothesis is not vacuous. The paper's standing $\alpha\ge3$ is replaced by each statement's own hypothesis $\alpha>0$.
• The minimum-norm minimizer and the Tikhonov points are predicates on a candidate point; where the curve $\epsilon\mapsto x_\epsilon$ is needed, a selection that minimizes $g+\frac\epsilon2\|\cdot\|^2$ for every $\epsilon>0$ is a hypothesis.
• $\int_{t_0}^{+\infty}\epsilon(t)/t\,dt=+\infty$ is stated as divergence of $T\mapsto\int_{t_0}^T\epsilon(t)/t\,dt$, never as an equation for a Bochner integral, which would take the value $0$ on a non-integrable function. Likewise, $\liminf\|x(t)-x^*\|=0$ is stated as "for every $\delta>0$, $\|x(t)-x^*\|<\delta$ for arbitrarily large $t$", not through Filter.liminf, whose value on an unbounded function is a default. $\sup<+\infty$ is boundedness from above of the image of $[t_0,+\infty)$, and $L^1$ membership is integrability on $[t_0,+\infty)$.
• Definitions needed: argmin, the minimum-norm minimizer, Tikhonov points, the General assumption, the solution predicate for (5), and $h_{x^*}$, $\dot h_{x^*}$. These duplicate objects of the other two missions of the series, which are drafted independently.
• Contributions welcome: proofs of the Tikhonov-curve facts in Mathlib generality, the energy estimate of Lemma 4.1, a continuous l'Hospital/Cesàro lemma for weighted averages, and the goal itself.

Selected references

• R.I. Boţ, E.R. Csetnek, S.C. László, Tikhonov regularization of a second order dynamical system with Hessian driven damping, Mathematical Programming, 2021, https://doi.org/10.1007/s10107-020-01528-8 (cited from arXiv:1911.12845v2). https://arxiv.org/abs/1911.12845
• H. Attouch, Z. Chbani, H. Riahi, Combining fast inertial dynamics for convex optimization with Tikhonov regularization, Journal of Mathematical Analysis and Applications 457(2), 1065–1094, 2018. https://arxiv.org/abs/1602.01973
• H. Attouch, J. Peypouquet, P. Redont, Fast convex optimization via inertial dynamics with Hessian driven damping, Journal of Differential Equations 261(10), 5734–5783, 2016. https://arxiv.org/abs/1601.07113
• H. Attouch, Z. Chbani, J. Peypouquet, P. Redont, Fast convergence of inertial dynamics and algorithms with asymptotic vanishing viscosity, Mathematical Programming 168, 123–175, 2018. https://arxiv.org/abs/1507.04782
• W. Su, S. Boyd, E.J. Candès, A differential equation for modeling Nesterov's accelerated gradient method: theory and insights, Journal of Machine Learning Research 17(153), 1–43, 2016. https://arxiv.org/abs/1503.01243

Motivation

Peaceman–Rachford splitting separates the minimization of a sum into two proximal computations. That matters when each summand has a tractable proximal map but the sum does not. Douglas–Rachford splitting is its half-relaxed instance, and related iterations appear in the analysis of ADMM. Davis and Yin study how regularity of the two summands changes the rate at which a relaxed Peaceman–Rachford iteration approaches a fixed point. Their Section 4 distinguishes two arrangements: one summand may carry both strong convexity and a Lipschitz gradient, or the two properties may be shared between summands. The mixed case of Theorem 4.3 is the target of this mission. The authors describe it as a new result in the preprint version used here. Davis and Yin, §4, pp. 13–16.

Setting

Let $H$ be a real Hilbert space and let $f,g:H\to(-\infty,+\infty]$ be proper, closed, convex functions. For a step size $\gamma>0$, the proximal map $P_f=\operatorname{prox}_{\gamma f}$ takes a point $z$ to the minimizer of $f(x)+\|x-z\|^2/(2\gamma)$; define $P_g$ similarly. The reflection associated with a proximal map $P$ is $R_P=2P-I$. The Peaceman–Rachford operator is $T_{\rm PRS}=R_{P_f}\circ R_{P_g}$. Its relaxed iteration, Algorithm 1 of the paper, is

$$z^{k+1}=(1-\lambda_k)z^k+\lambda_k T_{\rm PRS}(z^k), \qquad 0<\lambda_k\le1.$$

The analysis fixes a point $z^*$ satisfying $T_{\rm PRS}(z^*)=z^*$ and writes $x^*=P_g(z^*)$. It also uses the two proximal points $x_g^k=P_g(z^k)$ and $x_f^k=P_f(R_{P_g}(z^k))$. Their selected proximal subgradients are $\widetilde\nabla g(x_g^k)=(z^k-x_g^k)/\gamma$ and $\widetilde\nabla f(x_f^k)=(R_{P_g}(z^k)-x_f^k)/\gamma$. Lemma 1.1 states that they belong to the respective convex subdifferentials and relates them to the iteration. These are selected subgradients, rather than arbitrary members of a set-valued subdifferential. Davis and Yin, §1.2 and Lemma 1.1, pp. 4–7.

A function is $\mu$-strongly convex when its convexity inequality improves by the quadratic term $\mu t(1-t)\|x-y\|^2/2$. A $(1/\beta)$-Lipschitz gradient means the function is real valued, differentiable, and $\|\nabla h(x)-\nabla h(y)\|\le\|x-y\|/\beta$, with $\beta>0$. Section 1.10 also permits the parameters $\mu$ and $\beta$ to be zero when a property is unavailable. Its auxiliary term $S_f(x,y)$ is the maximum of a squared-point-distance term weighted by $\mu_f/2$ and a squared-selected-subgradient-distance term weighted by $\beta_f/2$. The one-step bound (2.1) controls the sum of the corresponding terms for $f$ and $g$. Davis and Yin, (1.13)–(1.14) and (2.1), pp. 7–8.

Formalization targets

Mixed regularity: Theorem 4.3

Suppose $\mu,\beta>0$. The goal covers both arrangements: $f$ is $\mu$-strongly convex while $g$ has a $(1/\beta)$-Lipschitz gradient, and the arrangement with $f$ and $g$ exchanged. In either case the target is

$$\|z^{k+1}-z^*\|\le C(\lambda_k)\|z^k-z^*\|, \qquad C(t)=\left(1-\frac{4t}{3}\min\left\{\gamma\mu,\frac{\beta}{\gamma},1-t\right\}\right)^{1/2}.$$

The factor is the paper's explicit factor, for $0\le t\le1$. It is strictly below one for $0<t<1$ with fixed positive regularity parameters; at $t=1$ it equals one. A uniform geometric rate therefore requires relaxation parameters to remain in a range where these factors have a uniform upper bound below one. Davis and Yin, Theorem 4.3, p. 16.

Companion rates

Theorems 4.1 and 4.2 give the explicit factors when $g$, or respectively $f$, supplies both properties. Proposition 4.1 translates a fixed-point contraction sequence into bounds on the proximal points, selected subgradients, fixed-point residual, and objective error. These are companion items in the proposal. The milestone list follows Lemma 1.1, the auxiliary bound (2.1), and the displayed estimates (4.6) and (4.7) from the proof of the mixed case. Davis and Yin, Proposition 4.1 and Theorems 4.1–4.3, pp. 13–16.

Significance

Theorem 4.3 identifies a rate available when neither summand alone carries both regularity properties. The explicit dependence on $\gamma$, $\mu$, $\beta$, and $\lambda_k$ lets a reader compare choices of step size and relaxation without hiding them in an unspecified constant. Proposition 4.1 explains what a contraction in the auxiliary fixed-point variable says about the quantities used to judge an optimization method: proximal point error, residual, and objective error. The distinction between a per-step bound and a uniform geometric rate is part of the mathematical result.

The inequalities are proved in the source paper; the remaining work here is to give their statements and eventual proofs a machine-checked account. The shared setting makes the proximal operator, selected subgradients, and regularity terms explicit, so the same formal definitions can be used across this mission's estimates. The two external platform definitions cited by this development already provide the closed proper convex class, proximal-map predicate, smooth convex class, and subdifferential. The results in this proposal are open Lean statements until solvers supply proofs.

Difficulty

Strong convexity and Lipschitz smoothness control different quantities in the mixed case. The decrease inequality directly bounds distances between proximal points and between gradients, whereas the desired conclusion concerns $z^k-z^*$. Treating $T_{\rm PRS}$ merely as a nonexpansive operator yields no factor below one. The displayed estimate (4.7) expresses the fixed-point distance through the quantities controlled by (4.6). The proof must also identify the proximal subgradient of the smooth summand with its ordinary gradient, rather than assume that identity as part of the goal. Davis and Yin, proof of Theorem 4.3, p. 16.

Formalization scope

The Lean carrier is a complete real inner-product space, including the zero space. The functions take values in EReal to represent $+\infty$ outside their effective domains. The published IsProperClosedConvex, IsProx, IsSmoothConvex, and subgrad definitions supply the common function and operator conventions. A proximal map is passed as a function together with a proof that it minimizes the defining objective; under the standing assumptions it is uniquely determined. Every result involving a proximal map assumes $\gamma>0$. Run statements require $0<\lambda_k\le1$, and the mixed theorem requires $\mu,\beta>0$. The goal has no hypothesis asserting any of its decrease estimates or a contraction; those are results to establish.

The contraction factors use Real.sqrt. For the mixed and $f$-regular factors, the radicands are positive throughout the relaxation range. For the $g$-regular factor, its nonnegativity follows from the compatibility of positive strong convexity and gradient Lipschitz parameters on a nonzero space; the zero space makes the contraction statement trivial even if that radicand is negative. Objective errors use real values only at proximal points where both summands are finite. The milestone (2.1) uses the paper's selected proximal subgradients. Proving the link between smooth gradients and those subgradients, and reusable proximal-map inequalities, are welcome contributions beyond the four named milestones.

Selected references

• Damek Davis and Wotao Yin, Faster convergence rates of relaxed Peaceman–Rachford and ADMM under regularity assumptions, arXiv:1407.5210v3, 2015; published in Mathematics of Operations Research 42(3), 2017. Preprint, journal DOI.

Motivation

The convex feasibility problem asks for a point in the intersection of closed convex sets. It underlies signal and image reconstruction, phase retrieval relaxations, and the constraint handling inside larger optimization methods. The workhorse algorithms are projection methods: the method of alternating projections (MAP), going back to von Neumann for subspaces, and the Douglas–Rachford and Peaceman–Rachford splitting schemes, which can be applied to the indicator, distance or squared-distance functions of the sets.

Without a regularity condition on how the sets meet, projection methods can converge arbitrarily slowly. Bounded linear regularity, an error bound saying that a point close to every set is close to their intersection, is the standard condition under which linear rates are obtained (Bauschke and Borwein, SIAM Review, 1996, doi:10.1137/S0036144593251710). Section 5 of D. Davis and W. Yin, Faster convergence rates of relaxed Peaceman–Rachford and ADMM under regularity assumptions (arXiv:1407.5210v3, Math. Oper. Res. 42(3), 2017) applies relaxed Peaceman–Rachford splitting (PRS) to the squared distance functions of two sets and proves that, under bounded linear regularity, the distance of the iterates to the intersection contracts by an explicit factor at every step, with step sizes and relaxation parameters allowed to vary between iterations. MAP is the special case of step size $1/2$ and relaxation $1$.

Setting

Let $\mathcal H$ be a real Hilbert space and $C_f, C_g \subseteq \mathcal H$ closed convex sets with $C_f\cap C_g\neq\emptyset$. For $C\subseteq\mathcal H$ the distance function is $d_C(x)=\inf_{y\in C}\|x-y\|$, and $P_C$ is the metric projection, the nearest point of $C$. The problem is modelled with

$$f(x)=d^2_{C_f}(x),\qquad g(x)=d^2_{C_g}(x).$$

For $\gamma>0$ the proximal map $\mathbf{prox}_{\gamma h}(x)$ is the minimiser of $h(y)+\frac1{2\gamma}\|y-x\|^2$, and the reflection is $\mathbf{refl}_{\gamma h}=2\,\mathbf{prox}_{\gamma h}-I$. With two step sizes, $T^{\gamma_f,\gamma_g}_{\mathrm{PRS}}=\mathbf{refl}_{\gamma_f f}\circ\mathbf{refl}_{\gamma_g g}$, and $(T)_\lambda=(1-\lambda)I+\lambda T$.

Given $z^0$, step sizes $\gamma_{f,k},\gamma_{g,k}>0$ and relaxation parameters $\lambda_k\in(0,1]$, iteration (5.1) is

$$x_g^k=\mathbf{prox}_{\gamma_{g,k}d^2_{C_g}}(z^k),\quad x_f^k=\mathbf{prox}_{\gamma_{f,k}d^2_{C_f}}(2x_g^k-z^k),\quad z^{k+1}=z^k+2\lambda_k(x_f^k-x_g^k).$$

The pair $\{C_f,C_g\}$ is boundedly linearly regular (Definition 5.1) if for every $\rho>0$ there is $\mu_\rho>0$ such that

$$d_{C_f\cap C_g}(x)\le\mu_\rho\max\{d_{C_f}(x),d_{C_g}(x)\}\qquad\text{for all }x\in B(0,\rho),$$

where $B(0,\rho)$ is the open ball. The rate function is

$$C(\gamma_f,\gamma_g,\lambda,\mu)=\left(1-\frac{4\lambda\min\{\gamma_g/(2\gamma_g+1)^2,\ \gamma_f/(2\gamma_f+1)^2\}}{\mu^2\max\{16\gamma_g^2/(2\gamma_g+1)^2,\ 1\}}\right)^{1/2}.$$

Formalization targets

Goal: Theorem 5.1

Suppose $\{C_f,C_g\}$ is boundedly linearly regular, $(z^k)$ is generated by (5.1), every $z^j$ lies in $B(0,\rho)$, and $d_{C_f\cap C_g}\le\mu_\rho\max\{d_{C_f},d_{C_g}\}$ on $B(0,\rho)$. Then for all $k\ge0$

$$d_{C_f\cap C_g}(z^{k+1})\le C(\gamma_{f,k},\gamma_{g,k},\lambda_k,\mu_\rho)\,d_{C_f\cap C_g}(z^k), \tag{5.3}$$

and if $\overline C=\sup_j C(\gamma_{f,j},\gamma_{g,j},\lambda_j,\mu_\rho)<1$, the iterates converge to a point $x\in C_f\cap C_g$ with

$$\|z^k-x\|\le 2\,d_{C_f\cap C_g}(z^0)\prod_{i=0}^{k-1}C(\gamma_{f,i},\gamma_{g,i},\lambda_i,\mu_\rho). \tag{5.4}$$

Milestones

The milestones follow the paper's own argument: Proposition 5.1 Parts 2 and 3 (the gradient of $d_C^2$ and the closed form $\mathbf{prox}_{\gamma d_C^2}=\frac1{2\gamma+1}I+\frac{2\gamma}{2\gamma+1}P_C$); Proposition 5.2, the fundamental inequality specialised to squared distances; Propositions C.1 and C.2 (the fixed-point set of $T^{\gamma_f,\gamma_g}_{\mathrm{PRS}}$ is $C_f\cap C_g$, and the iterates are Fejér monotone for $\lambda_k\in(0,1]$); and the displays (5.5)–(5.7), (5.8), (5.10) and (5.11) of the proof of Theorem 5.1, ending with the one-step contraction at a single point.

Companion: Corollary 5.1

With $\gamma_{f,k}\equiv\gamma_{g,k}\equiv\frac12$ and $\lambda_k\equiv1$, iteration (5.1) is MAP, $z^{k+1}=P_{C_f}P_{C_g}z^k$, and under the assumptions of Theorem 5.1 it contracts the distance to $C_f\cap C_g$ by $(1-1/\mu_\rho^2)^{1/2}$ from the second iterate on, converging linearly to a point of the intersection.

Significance

Theorem 5.1 gives an explicit, iteration-wise linear rate for relaxed PRS on the two-set feasibility problem with step sizes and relaxation parameters that may change from step to step. Its specialisation, Corollary 5.1, recovers linear convergence of MAP under bounded linear regularity with rate $(1-1/\mu_\rho^2)^{1/2}$, better than the $(1-1/(8\mu^2))^{1/2}$ derived in earlier work on cyclic projections. Remark 5.1 of the paper notes that Douglas–Rachford on indicator functions is a limiting case of the scheme as the step sizes grow, though the rate degenerates in that limit.

The results are proved in the paper; the proof of the convergence statement cites Bauschke–Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Spaces, Theorem 5.12 (Fejér monotone sequences whose distance to a closed convex set decays converge strongly). No machine-checked proof of these statements was found on Prove2Me (search of 2026-10-07). A formalization produces a checked chain from the closed form of the prox of $d_C^2$ through the fundamental inequality to the linear rate, and with it reusable facts about squared distance functions, metric projections onto closed convex sets in Hilbert space, and Fejér monotone sequences.

Difficulty

The obvious argument would contract the distance to $C_f\cap C_g$ directly, but one PRS step measures only $d_{C_g}$ at $z$ and $d_{C_f}$ at the reflected point $\mathbf{refl}_{\gamma_g g}(z)$, never $d_{C_f}(z)$. Transferring the second distance back to $z$ costs a factor $\max\{c_1^2,1\}$ with $c_1=4\gamma_g/(2\gamma_g+1)$, and only after that does the regularity inequality turn $\max\{d_{C_f}(z),d_{C_g}(z)\}$ into $d_{C_f\cap C_g}(z)$. The fundamental inequality (5.2) itself rests on the prox-gradient calculus of the earlier sections, specialised to functions whose gradients $2(I-P_C)$ are Lipschitz. Passing from a contracting distance to convergence of the iterates to one point needs Fejér monotonicity and a completeness argument; a contracting distance alone does not give a limit.

Formalization scope

Everything lives in a real Hilbert space (InnerProductSpace ℝ H with CompleteSpace H); $d_C$ is Mathlib's Metric.infDist. The squared distance is viewed as an extended-real function so that the published proximal predicate ThreeOpSplitting.ConvexRates.IsProx applies, and prox maps are arbitrary maps satisfying it (unique here). Metric projections are maps satisfying the published RandomGradFree.Nonsmooth.IsMetricProjection. Step sizes and prox maps are indexed by the iteration. The hypotheses of the goal are: $C_f,C_g$ closed and convex with nonempty intersection; bounded linear regularity of the pair; $\gamma_{f,k},\gamma_{g,k}>0$; $\lambda_k\in(0,1]$; iterates in an open ball $B(0,\rho)$ on which the regularity inequality holds with $\mu_\rho>0$. The range $\lambda_k\in(0,1]$ is Algorithm 1's standing assumption, used in the proof but not restated in (5.1). The condition $\overline C<1$ is stated as a uniform bound $c<1$ on all rate factors. The product in (5.4) runs over $i=0,\dots,k-1$; the page prints the upper index $k$, which claims an extra factor that the proof does not give. In the proof's (5.9)–(5.10) the page swaps $\gamma_f$ and $\gamma_g$ relative to (5.2); the milestone states the pairing that (5.2) produces. Corollary 5.1's improved rate is stated for $k\ge1$, since it uses $z^k\in C_f$.

A trivializing formalization is ruled out: the goal mentions none of the proof's intermediate quantities, takes the prox maps as genuine minimisers of the squared-distance objective rather than chosen points, and measures distances to a nonempty intersection, so the junk value of infDist on the empty set never arises.

A complete development needs the closed form of $\mathbf{prox}_{\gamma d^2_C}$, differentiability of $d_C^2$, the fundamental inequality of the paper's §1 specialised to these functions, and the convergence theorem for Fejér monotone sequences. The last three are reusable well beyond this mission. Proofs of any milestone, and of Mathlib-level lemmas on metric projections, are welcome.

Selected references

• D. Davis, W. Yin, Faster convergence rates of relaxed Peaceman–Rachford and ADMM under regularity assumptions, Mathematics of Operations Research 42(3), 2017; arXiv version used here: arXiv:1407.5210v3.
• H. H. Bauschke, J. M. Borwein, On projection algorithms for solving convex feasibility problems, SIAM Review 38(3), 1996. doi:10.1137/S0036144593251710
• H. H. Bauschke, P. L. Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Spaces, Springer, 2011. doi:10.1007/978-1-4419-9467-7
• P.-L. Lions, B. Mercier, Splitting algorithms for the sum of two nonlinear operators, SIAM Journal on Numerical Analysis 16(6), 1979. doi:10.1137/0716071