Numerical Analysis

Motivation

Most equations that arise in applications cannot be solved in closed form: the age of the Moon from a radioactive-decay balance, the deflection of a clamped beam, the equilibrium of a catenary, all reduce to solving $f(x) = 0$ for a function $f$ with no algebraic inverse. A first course in numerical methods therefore opens with root finding: constructive procedures that produce a sequence of approximations $x_0, x_1, x_2, \\dots$ together with a theorem saying that the sequence converges to a root and how fast.

This mission is the first in a series formalizing the lecture notes Métodos Numéricos by Sergio Roberto de Freitas (Departamento de Computação e Estatística, UFMS, 2000). It covers Chapter 3, Zeros de Funções, whose capstone is the local convergence of Newton's method.

Setting

Let $f : \\mathbb{R} \\to \\mathbb{R}$ and let $\\bar{x}$ be a zero of $f$, i.e. $f(\\bar{x}) = 0$. The chapter studies three constructions.

Bisection. Starting from an interval $(a_0, b_0)$ with $f(a_0) < 0 < f(b_0)$, put $x_{k+1} = (a_k + b_k)/2$ and keep the half of $(a_k, b_k)$ on whose endpoints $f$ still changes sign. The width of the bracketing interval is halved at every step.

Linear iteration (MIL). Rewrite $f(x) = 0$ as a fixed-point equation $x = g(x)$ and iterate $x_{n+1} = g(x_n)$ from an arbitrary $x_0$.

Newton's method. Take the particular iteration function

$$g(x) = x - \\frac{f(x)}{f'(x)}, \\qquad x_{n+1} = x_n - \\frac{f(x_n)}{f'(x_n)},$$

obtained by truncating the Taylor expansion of $f$ at $x_n$ after the linear term.

A sequence $x_n \\to \\alpha$ has order of convergence $p$ when $|e_{n+1}| / |e_n|^p$ tends to a finite constant, where $e_n = x_n - \\alpha$; $p = 1$ is linear and $p = 2$ quadratic convergence.

Target

The goal theorem is the local convergence of Newton's method at a simple zero. If $f$ is twice differentiable on an open interval $(a, b)$ containing $\\bar{x}$, its second derivative is continuous there, and $f'$ never vanishes on $(a, b)$, then there is $h > 0$ such that

$$x_0 \\in [\\bar{x} - h, \\bar{x} + h] \\implies x_n \\to \\bar{x},$$

where $x_n$ is the Newton sequence started at $x_0$. This is the formal content of the chapter's statement that Newton's method converges provided the initial guess is chosen close enough to the root.

The milestones are the supporting results of the chapter: the error bound and convergence of bisection, the convergence of the linear iterative method under a derivative bound $|g'| \\le L < 1$, the a posteriori estimate $|\\bar{x} - x_n| \\le \\frac{L}{1-L}|x_n - x_{n-1}|$, and the quadratic order of Newton's method at a simple zero.

Significance

Bisection, fixed-point iteration and Newton's method are the three root finders every numerical-analysis course starts with, and the three convergence theorems above are what justifies using them. The a posteriori estimate is what turns the iteration into an algorithm with a stopping criterion: it bounds the distance to the root by a quantity the program can measure. The quadratic order statement explains the observed doubling of correct digits per step at a simple zero and its loss at a multiple zero.

On the formalization side, Mathlib already has a fixed-point theorem for contractions on complete spaces and the mean value theorem, but not the statements in the form used in numerical analysis: bisection with its explicit $2^{-n}$ bracket, the $\\frac{L}{1-L}$ a posteriori bound, or Newton's local convergence and quadratic rate stated for the concrete iteration sequence. This mission asks for those.

Difficulty

The delicate point in all three theorems is that the iterates must be known to stay in the region where the hypotheses hold; the informal proofs assume this silently. For the linear iterative method the mission therefore states the invariance hypothesis explicitly ($g$ maps the closed interval into itself). For Newton's method, no such hypothesis is given: the existence of a neighbourhood of $\\bar{x}$ that the iteration preserves is part of what must be proved, and it comes from the continuity of $g'$ together with $g'(\\bar{x}) = 0$. The quadratic-order milestone also has to handle the degenerate possibility $x_n = \\bar{x}$, which would put a zero in the denominator; it is excluded by hypothesis.

Formalization scope

Everything is over the real numbers. The iterations are given as explicit recursive sequences: the bisection construction returns the bracketing pair $(a_n, b_n)$ and its midpoint, and there are separate sequences for the fixed-point and Newton iterations. Newton's method takes the derivative as a separate function argument $f'$, tied to $f$ by a hypothesis of the form "$f$ has derivative $f'(x)$ at every $x$ of the interval"; this avoids relying on any junk value of a derivative operator where $f$ fails to be differentiable. Division is total, so a step at a point where $f'$ vanishes would leave the iterate unchanged; the hypotheses exclude this inside the interval.

Two deliberate deviations from the source are worth flagging for the auditor. First, Proposição 3.5.1 assumes $f'(x) \\neq 0$ only for $x \\neq \\bar{x}$, but its proof divides by $f'(\\bar{x})^2$; the formal statement assumes $f' \\neq 0$ on the whole interval, so $\\bar{x}$ is a simple zero. Second, the convergence statement for the linear iterative method adds the hypothesis that $g$ maps the closed interval into itself, without which the iterates may leave the region where the derivative bound is assumed.

Selected references

• S. R. Freitas, Métodos Numéricos, Departamento de Computação e Estatística, Centro de Ciências Exatas e Tecnologia, Universidade Federal de Mato Grosso do Sul, 2000. Chapter 3, Zeros de Funções, pp. 33–70. (Course notes supplied with this mission.)

Motivation

A table of values is all one has of many functions: measurements, tabulated physical constants, the output of an expensive simulation. Interpolation reconstructs a function between tabulated points by passing a polynomial through them, and it underlies much of the rest of numerical analysis — quadrature rules, finite-difference formulas and predictor-corrector schemes for differential equations are all obtained by integrating or differentiating an interpolating polynomial. What makes the reconstruction trustworthy is a formula for the error committed away from the nodes.

This mission is the fifth in a series formalizing the lecture notes Métodos Numéricos by Sergio Roberto de Freitas (UFMS, 2000). It covers Chapter 7, Interpolação.

Setting

Let $x_0 < x_1 < \\dots < x_n$ be distinct nodes and $f_i = f(x_i)$ the tabulated values of a function $f$. The Lagrange form of the interpolating polynomial is

$$P_n(t) = \\sum_{i=0}^{n} f_i L_i(t), \\qquad L_i(t) = \\prod_{j \\neq i}\\frac{t - x_j}{x_i - x_j},$$

so that $P_n(x_i) = f_i$ for every $i$, and $P_n$ is the unique polynomial of degree at most $n$ with this property.

The interpolation error at a point $x$ that is not a node is the difference $f(x) - P_n(x)$.

Target

The goal theorem is Proposição 7.3.1: if $f$ is $n+1$ times differentiable on $(a,b)$ and all nodes lie in $(a,b)$, then for every $x \\in (a,b)$ different from all nodes there exists $\\xi \\in (a,b)$ with

$$f(x) = P_n(x) + \\left(\\prod_{k=0}^{n}(x - x_k)\\right)\\frac{f^{(n+1)}(\\xi)}{(n+1)!}.$$

The milestones are the other results of the chapter needed to make sense of it: the existence and uniqueness of the interpolating polynomial of degree at most $n$ through $n+1$ points with distinct abscissas (Proposição 7.2.1), the fact that the Lagrange formula does interpolate the data, and the practical error bound of Proposição 7.3.2, $|f(x) - P_n(x)| \\le \\frac{M}{(n+1)!}\\prod_k |x - x_k|$ when $|f^{(n+1)}| \\le M$.

Significance

The error formula is the reason interpolation is a numerical method and not just a curve-drawing device: it shows that the error is governed by two independent factors, the smoothness of $f$ through $f^{(n+1)}$ and the geometry of the nodes through the product $\\prod_k(x - x_k)$. Everything downstream follows from it — the $h^2$ and $h^4$ error terms of the trapezoidal and Simpson rules in the next mission are obtained by integrating exactly this expression, and the choice of Chebyshev nodes is the attempt to make the product factor small.

Difficulty

The standard proof introduces the auxiliary function $F(t) = f(t) - P_n(t) - K\\prod_k(t - x_k)$ with $K$ chosen so that $F(x) = 0$, and then applies Rolle's theorem $n+1$ times to conclude that $F^{(n+1)}$ vanishes somewhere. Formalizing the repeated application of Rolle's theorem, keeping track of the $n+2$ distinct zeros and the nested intervals they generate, is the substance of the work; it is an induction that has to be organized carefully rather than a computation.

Formalization scope

The nodes are given as a strictly increasing family $x_0 < \\dots < x_n$ of $n+1$ reals lying in the open interval $(a,b)$, and the interpolating polynomial is the explicit Lagrange sum rather than an abstract polynomial: at a node it is defined by the same formula, whose factors then include $0/0$ contributions unless the nodes are distinct, which is why the distinctness hypothesis appears in the interpolation milestone. Smoothness is expressed as continuous differentiability of order $n+1$ on the open interval $(a,b)$, and the derivative appearing in the error term is the iterated derivative computed within that set; since the set is open this agrees with the ordinary $(n+1)$-st derivative. The evaluation point $x$ is assumed to lie in $(a,b)$ and to differ from every node; the point $\\xi$ is asserted to exist in $(a,b)$, with no claim of uniqueness or of any relation to $x$ beyond membership in the interval. The uniqueness milestone is stated with Mathlib's polynomial type and the degree bound $\\deg \\le n$, which includes the zero polynomial.

Selected references

• S. R. Freitas, Métodos Numéricos, Departamento de Computação e Estatística, Universidade Federal de Mato Grosso do Sul, 2000. Chapter 7, Interpolação, pp. 133–161. (Course notes supplied with this mission.)

Why the image function matters for ADMM

The alternating direction method of multipliers (ADMM) solves problems of the form $\min_{x,z} f(x)+g(z)$ subject to $Ax+Bz=b$. Themelis and Patrinos (arXiv:1709.05747v4, SIAM J. Optim. 30(1):149–181, 2020) show that, for nonconvex $f$ and $g$, ADMM is the Douglas–Rachford splitting method applied to a reformulated problem whose second term is the image function $\varphi_2=(Bg)$, $(Bg)(s)=\inf\{g(z)\mid Bz=s\}$. Their convergence theory for DRS requires the nonsmooth term to be proper and lower semicontinuous (Assumption IIa4), so whether $(Bg)$ inherits lower semicontinuity from $g$ decides whether the theory applies to a given ADMM instance.

Image functions, also called epi-compositions or infimal post-compositions, are classical objects of variational analysis (Rockafellar and Wets, Variational Analysis; Bauschke and Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Spaces; Auslender and Teboulle, Asymptotic Cones and Functions). For convex $g$ the question of when $(Bg)$ is closed has a long history through recession functions; Section 5.4.1 of the paper gives a characterization that needs no convexity.

Setting

Let $\overline{\mathbb R}=\mathbb R\cup\{+\infty\}$. For $g:\mathbb R^n\to\overline{\mathbb R}$ the domain is $\operatorname{dom}g=\{z\mid g(z)<\infty\}$; $g$ is lower semicontinuous (lsc) if its epigraph is closed, and level bounded if every set $\{z\mid g(z)\le\alpha\}$, $\alpha\in\mathbb R$, is bounded.

Fix a matrix $B\in\mathbb R^{p\times n}$. The image function of $g$ under $B$ (Definition 5.1) is

$$(Bg)(s)=\inf_{z\in\mathbb R^n}\{g(z)\mid Bz=s\},\qquad s\in\mathbb R^p,$$

with values in $[-\infty,+\infty]$. Its domain is $B\operatorname{dom}g$: outside that set the infimum is over values $+\infty$ only.

Condition (5.5) at a point $\bar z\in\operatorname{dom}g$ compares the best value of $g$ on the affine fibre $\bar z+\ker B$ with the values of $g$ far away along directions that $B$ nearly annihilates:

$$\liminf_{\substack{\|d\|\to\infty\\ Bd\to0}} g(\bar z+d)\;\ge\;\inf_{d\in\ker B} g(\bar z+d).$$

The right side equals $(Bg)(B\bar z)$.

Formalization targets

Goal: Theorem 5.11, on $B\operatorname{dom}g$

For every lsc $g:\mathbb R^n\to\overline{\mathbb R}$ and every $B\in\mathbb R^{p\times n}$,

$$(Bg)\ \text{is lsc at every } \bar s\in B\operatorname{dom}g\quad\Longleftrightarrow\quad (5.5)\ \text{holds at every } \bar z\in\operatorname{dom}g.$$

Milestones

1. The right side of (5.5) is $(Bg)(B\bar z)$.
2. The "if" direction: (5.5) on $\operatorname{dom}g$ gives lower semicontinuity of $(Bg)$ at every point of $B\operatorname{dom}g$.
3. The "only if" direction: lower semicontinuity of $(Bg)$ on $B\operatorname{dom}g$ forces (5.5) on $\operatorname{dom}g$.

Companion results

• The "In particular" clause of Theorem 5.11, as printed: if $g$ is lsc and level bounded, then $(Bg)$ is lsc on all of $\mathbb R^p$.
• The properness clause of Proposition 5.10: if $g$ is proper and lsc and, for some $\beta>0$, $\min_z\{g(z)+\frac\beta2\|Bz-s\|^2\}$ is attained for every $s$, then $(Bg)$ is proper.

Significance

Theorem 5.11 turns Assumption IIa4, a property of a derived function, into a checkable condition on $g$ and $B$. Its level-bounded corollary covers every lsc $g$ with bounded sublevel sets, and the characterization explains why the pathological example (5.4) of the paper, in which $g$ is even $C^1$, fails: a minimizing sequence on the fibre escapes to infinity while $B$ applied to it converges. Combined with the convergence theorem for ADMM in the same paper (Theorem 5.6), this identifies the ADMM instances to which the tight DRS stepsize theory applies.

The theorem is proved in the paper. As far as is known, no machine-checked proof of it, or of any closedness criterion for image functions, exists in Lean or in Mathlib. A formalization fixes the exact range of the equivalence, which the printed statement gets wrong (see Formalization scope), and supplies reusable infrastructure: the image function in EReal, and liminfs along the filter "$\|d\|\to\infty$, $Bd\to0$".

Difficulty

Lower semicontinuity of $(Bg)$ at $\bar s$ asks that values of $g$ on fibres $\{Bz=s\}$, $s\to\bar s$, not drop below $(Bg)(\bar s)$. The natural argument picks near-minimizers $z_k$ on these fibres and passes to a limit; this works only when the $z_k$ stay bounded, and nothing in lsc of $g$ alone prevents them from escaping to infinity. The difficulty is to account exactly for the unbounded case, which must be compared with the fibre through $\bar s$ by a translation $z_k-\bar z$, and that requires a point $\bar z\in\operatorname{dom}g$ with $B\bar z=\bar s$. The converse requires turning a failure of (5.5), stated as a filter liminf, into a sequence witnessing the failure of lower semicontinuity.

Formalization scope

• $\mathbb R^n$ and $\mathbb R^p$ are EuclideanSpace ℝ (Fin n) and EuclideanSpace ℝ (Fin p); $B$ is a continuous linear map between them. Every dimension, including $0$, is allowed.
• $\overline{\mathbb R}$ is EReal, and $g$ is assumed never to take the value $-\infty$. The image function is the EReal infimum $\inf_{z,\ Bz=s} g(z)$, which may be $-\infty$ and is $+\infty$ off $B\operatorname{dom}g$. Properness of $g$ is not assumed in the goal: if $\operatorname{dom}g=\emptyset$ both sides hold.
• Lower semicontinuity at a point is Mathlib's LowerSemicontinuousAt, the neighbourhood form of the page's sequential level-set criterion.
• The liminf in (5.5) is taken in the complete lattice $[-\infty,+\infty]$ along the meet of the cobounded filter ($\|d\|\to\infty$) and the preimage under $B$ of the neighbourhood filter of $0$. When no such $d$ exist (e.g. $B$ injective) the filter is trivial, the liminf is $+\infty$, and (5.5) holds, as an empty liminf does on the page.
• Correction. As printed, Theorem 5.11 asserts lsc of $(Bg)$ everywhere, and its "if" direction is then false. Take $n=2$, $p=1$, $B=[1\ 0]$, and $g$ the indicator of the closed set $\{(x,y):x>0,\ xy\ge1\}$. Condition (5.5) holds at every $\bar z\in\operatorname{dom}g$, since its right side is $0$ and $g\ge0$; but $(Bg)(s)=0$ for $s>0$ and $+\infty$ for $s\le0$, so $(Bg)$ is not lsc at $0$. The proof begins "given $\bar s\in\operatorname{dom}(Bg)$" and covers exactly $B\operatorname{dom}g$; both directions establish the corrected equivalence, which is the goal. The "In particular" clause is true as printed and is stated with lsc everywhere. The lsc clause of Proposition 5.10 fails for the same reason and is not formalized.
• Assumption IIa2 enters Proposition 5.10 only through the solvability of $\min_z\{g(z)+\frac\beta2\|Bz-s\|^2\}$ for every $s$.
• A trivializing formalization is ruled out: the goal assumes only lsc of $g$ and linearity of $B$, as on the page, and is not restricted to proper or level-bounded $g$ (that is the separate companion).
• Contributions welcome: proofs of the milestones, a general Mathlib-style lemma that the infimum of an lsc function over the fibres of a linear map is lsc on its domain under (5.5), and the closedness of image functions under level-boundedness.

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. Preprint arXiv:1709.05747v4. https://arxiv.org/abs/1709.05747, https://doi.org/10.1137/18M1163993
• R. T. Rockafellar, R. J.-B. Wets, Variational Analysis, Grundlehren 317, Springer, 1998 (corrected printings 2009, 2011). https://doi.org/10.1007/978-3-642-02431-3
• H. H. Bauschke, P. L. Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Spaces, 2nd ed., Springer, 2017. https://doi.org/10.1007/978-3-319-48311-5
• A. Auslender, M. Teboulle, Asymptotic Cones and Functions in Optimization and Variational Inequalities, Springer, 2003 (cited as 2002 in the paper). https://doi.org/10.1007/b97594

Motivation

The alternating direction method of multipliers (ADMM) solves linearly constrained problems

$$\underset{(x,z)\in\mathbb R^m\times\mathbb R^n}{\text{minimize}}\ f(x)+g(z)\quad\text{subject to}\quad Ax+Bz=b,\tag{1.2}$$

by alternating a minimization in $x$, a minimization in $z$ and a multiplier update. It is one of the standard methods for distributed optimization, statistical learning and signal processing, and in practice it is often applied to nonconvex $f$ or $g$ (see Boyd et al., Distributed optimization and statistical learning via the alternating direction method of multipliers, 2011, doi:10.1561/2200000016).

For convex problems, ADMM with $\lambda=1$ is classically Douglas–Rachford splitting (DRS) applied to the Fenchel dual of (1.2) (Gabay 1983). The dual reading needs convexity. Themelis and Patrinos (arXiv:1709.05747v4, §5; SIAM J. Optim. 30(1), 2020, doi:10.1137/18M1163993) prove a primal equivalence instead: for arbitrary $f,g$ and for every relaxation parameter $\lambda$, one ADMM update is one DRS step on a reformulation of (1.2) built from image functions. A primal equivalence for $\lambda=1$ was noted earlier by Bauschke and Koch (2015, Rem. 3.14) for $A=-B=I$, and Yan and Yin (2016, Thm. 1) showed that this case loses no generality. This mission formalizes the general statement and the convergence theorem for nonconvex ADMM that the paper derives from it.

Setting

All spaces are Euclidean, $\overline{\mathbb R}=\mathbb R\cup\{+\infty\}$, $A\in\mathbb R^{p\times m}$, $B\in\mathbb R^{p\times n}$, $b\in\mathbb R^p$. A function is proper if it never takes the value $-\infty$ and is finite somewhere. The proximal mapping of $h$ with stepsize $\gamma>0$ is the set

$$\operatorname{prox}_{\gamma h}(s)=\operatorname*{arg\,min}_{w}\Big\{h(w)+\tfrac1{2\gamma}\|w-s\|^2\Big\}.$$

For $h:\mathbb R^n\to\overline{\mathbb R}$ and a linear map $C$, the image function is

$$(Ch)(s)=\inf_{x}\{h(x)\mid Cx=s\}\in[-\infty,+\infty].$$

The $\beta$-augmented Lagrangian of (1.2) is $\mathcal L_\beta(x,z,y)=f(x)+g(z)+\langle y,Ax+Bz-b\rangle+\frac\beta2\|Ax+Bz-b\|^2$. One ADMM update with penalty $\beta>0$ and relaxation $\lambda>0$ maps $(x,y,z)$ to $(x^+,y^+,z^+)$:

$$y^{+/2}=y-\beta(1-\lambda)(Ax+Bz-b),\quad x^+\in\operatorname*{arg\,min}\mathcal L_\beta(\cdot,z,y^{+/2}),\quad y^+=y^{+/2}+\beta(Ax^++Bz-b),\quad z^+\in\operatorname*{arg\,min}\mathcal L_\beta(x^+,\cdot,y^+).$$

One DRS step with stepsize $\gamma$ and relaxation $\lambda$ on $\varphi_1+\varphi_2$ maps $s$ to $s^+=s+\lambda(v-u)$, where $u\in\operatorname{prox}_{\gamma\varphi_1}(s)$ and $v\in\operatorname{prox}_{\gamma\varphi_2}(2u-s)$. The reformulation of (1.2) is

$$\underset{s\in\mathbb R^p}{\text{minimize}}\ (Af)(s)+(Bg)(b-s),\tag{5.1}$$

that is, $\varphi_1=(Af)$ and $\varphi_2=(Bg)(b-\cdot)$. The regular subdifferential is $\hat\partial$, the limiting one $\partial$.

Formalization targets

Goal: primal equivalence (Theorem 5.5)

Let $f,g$ be proper, $\beta>0$ large enough that every ADMM subproblem has a solution, and $\lambda>0$. With $s=Ax-y/\beta$, $u=Ax$, $v=b-Bz$ and $s^+,u^+,v^+$ the same expressions at $(x^+,y^+,z^+)$, and $\gamma=1/\beta$,

$$s^+=s+\lambda(v-u),\qquad u^+\in\operatorname{prox}_{\gamma\varphi_1}(s^+),\qquad v^+\in\operatorname{prox}_{\gamma\varphi_2}(2u^+-s^+),$$

together with $\varphi_1(u^+)=f(x^+)$, $\varphi_2(v^+)=g(z^+)$, $-y^+\in\hat\partial\varphi_1(u^+)$, $-A^\top y^+\in\hat\partial f(x^+)$ and $\operatorname{dist}(-B^\top y^+,\hat\partial g(z^+))\le\beta\|B\|\|Ax^++Bz^+-b\|$.

Milestones

Proposition 5.2(i)–(iii) (properness of $(Ch)$, $(Ch)(Cx_\beta)=h(x_\beta)$, $\operatorname{prox}_{(Ch)/\beta}\supseteq CX_\beta$), Proposition 5.3 ($C^\top\hat\partial(Ch)(\bar s)\subseteq\hat\partial h(\bar x)$), the inclusions (5.3a) and (5.3b), and the identity $s+\lambda(v-u)=s^+$.

Companions

Theorem 5.5(vi): if $A$ is surjective, $(Af)$ is $L$-smooth and $\beta>L$, then $\varphi^{\mathrm{DR}}_\gamma(s^+)=\mathcal L_\beta(x^+,z^+,y^+)$. Theorem 5.6: under Assumption II and $(\gamma,\lambda)$ in the ranges of Theorem 4.1 with $\beta=1/\gamma$, $\mathcal L_\beta$ decreases by $\frac{c\lambda^2}{(1+\gamma L)^2}\|Ax^k+Bz^k-b\|^2$ per step, $\min_{i\le k}\|Ax^i+Bz^i-b\|=o(1/\sqrt k)$, cluster points are KKT points with common cost, and the iterates are bounded when $\Phi$ is level bounded.

Significance

Theorem 5.5 reduces nonconvex ADMM to nonconvex DRS: any descent, convergence or rate statement for DRS on $\varphi_1+\varphi_2$ becomes one for ADMM, with the Douglas–Rachford envelope becoming the augmented Lagrangian (5.5(vi)). This is how the paper obtains Theorem 5.6, with explicit penalty bounds for every $\lambda\in(0,4)$ permitted by Theorem 4.1, under assumptions weaker than earlier nonconvex ADMM analyses (Li and Pong, 2015, doi:10.1137/140998135, require $A=I$ and coercivity).

The results are proved in the paper. No machine-checked proof of a primal ADMM–DRS equivalence is known; the platform holds formalizations of proximal ADMM variants and of convex DRS for monotone operators, but not of image functions. The mission produces a reusable library of image functions and their proximal and subdifferential calculus, and a precise statement of relaxed ADMM in the nonconvex, extended-real setting.

Difficulty

The individual steps are short, but the objects are extended-real and set-valued. An image function may equal $-\infty$, and its infimum need not be attained; that $(Af)$ is proper and that the infimum is attained at $x^+$ must be derived from the solvability of the ADMM subproblems, not assumed. The proximal mappings are multivalued, so the equivalence is an inclusion between argmin sets, and the algebra relating $y^{+/2}$, $y^+$ and $s^+$ has to be matched to the relaxation exactly. The subgradient claims require regular subdifferentials of extended-real functions, and (v) requires the optimality condition of the $z$-subproblem, whose quadratic term couples $z^+$ with the residual. Theorem 5.6 additionally rests on the whole DRS convergence analysis of §4.

Formalization scope

Points are EuclideanSpace ℝ (Fin _), matrices are continuous linear maps, $A^\top$ is ContinuousLinearMap.adjoint, and $\|B\|$ is the operator norm. Extended-real functions are EReal-valued; the image function imageFn is the EReal infimum over $\{x\mid Cx=s\}$, so it is $+\infty$ off the range and may be $-\infty$. Proximal mappings and $X_\beta$ are sets, and an ADMM update (IsADMMStep) is any selection of minimizers in the order $y^{+/2},x^+,y^+,z^+$. The two subproblems minimize the terms of $\mathcal L_\beta$ that depend on the minimized variable, because the dropped constant $g(z)$ may be $+\infty$ at a starting point and would make every $x$ a minimizer. "Large enough penalty" is SubproblemsSolvable. The regular and limiting subdifferentials are the published IsRegularSubgrad and LimitingSubdiff.

Added hypotheses: $f$ and $g$ proper in Theorem 5.5 and 5.5(vi), and $h$ proper in Proposition 5.2(i). The page omits them, and they are necessary: for $h\equiv+\infty$ every $X_\beta(s)$ is the whole space while $(Ch)\equiv+\infty$ is not proper. Claim (v) is stated as the existence of $w\in\hat\partial g(z^+)$ with $\|-B^\top y^+-w\|\le\beta\|B\|\|Ax^++Bz^+-b\|$, which is stronger than the distance bound and avoids the distance to an empty set. Theorem 5.6 uses the constants $c$ of Theorem 4.1 as derived in its proof ((4.8), (4.10)); the printed (4.3) has the sign of one entry of a max reversed and the printed (4.5) an extra factor $1/\lambda$, and both are then false (counterexamples in the companion mission on DRS convergence). Theorem 5.6 adds $L>0$; its run is indexed so that the paper's $(x^k,y^k,z^k)$ is the term $k+1$, with the starting triple at index $0$.

The goal does not assume the DRS relations or (5.3a)–(5.3b), and $s^+,u^+,v^+$ are defined from the ADMM iterates by (5.2), not chosen independently; a formalization that quantifies over an arbitrary prox selection, or that hard-codes $\lambda=1$, is a different statement.

Needed infrastructure: image functions and their proximal calculus, regular subgradients of extended-real functions under linear maps, and, for Theorem 5.6, the DRS convergence theory of §4 (formalized in the companion mission I of this series). Proofs of individual milestones are welcome independently.

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
• H. H. Bauschke, V. R. Koch, Projection methods: Swiss army knives for solving feasibility and best approximation problems with halfspaces, in Infinite Products of Operators and Their Applications, Contemp. Math. 636, pp. 1–40, AMS, 2015.
• M. Yan, W. Yin, Self equivalence of the alternating direction method of multipliers, in Splitting Methods in Communication, Imaging, Science, and Engineering, pp. 165–194, Springer, 2016.
• G. Li, T. K. Pong, Global convergence of splitting methods for nonconvex composite optimization, SIAM J. Optim. 25(4):2434–2460, 2015. doi:10.1137/140998135
• 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), 2011. doi:10.1561/2200000016
• D. Gabay, Applications of the method of multipliers to variational inequalities, in Augmented Lagrangian Methods, Studies in Mathematics and Its Applications 15, pp. 299–331, Elsevier, 1983.

Motivation

Accelerated first-order methods, beginning with Nesterov's 1983 scheme, reach the optimal iteration complexity among gradient methods for smooth convex minimization, but their classical derivations ("estimate sequences") give little insight into why acceleration works. A line of work started by Su, Boyd and Candès (arXiv:1503.01243) reads accelerated methods as discretizations of second-order ordinary differential equations and proves their rates with Lyapunov functions: a quantity that decreases along the continuous trajectory and, after discretization, along the iterates.

Luo and Chen (arXiv:1909.03145, Math. Program. 2022) introduce the NAG flow, a first-order system in two variables $(x, v)$ with a time-varying parameter $\gamma$, and derive several accelerated methods from it as numerical ODE solvers. Their §4 treats the simplest discretization, the fully implicit (backward Euler) scheme. It is not a practical method by itself, since each step solves a nonlinear equation, but it is the case in which the continuous Lyapunov argument transfers to the discrete level almost word for word, and it is the benchmark for the semi-implicit and explicit schemes of §§5–7. This mission formalizes its convergence theorem, Theorem 4.1.

Setting

Let $V$ be a real Hilbert space with inner product $(\cdot,\cdot)$ and norm $\|\cdot\|$; the duality pairing $\langle\cdot,\cdot\rangle$ is identified with the inner product. Let $f:V\to\mathbb R$ and $\mu\ge 0$. The function $f$ belongs to the class $\mathcal S^1_\mu$ if it is continuously differentiable with gradient $\nabla f$ and satisfies the $\mu$-convexity inequality

$$f(x)-f(y)-\langle\nabla f(y),x-y\rangle\ \ge\ \frac\mu2\|x-y\|^2\qquad\forall x,y\in V .$$

For $\mu=0$ this is convexity, for $\mu>0$ strong convexity. A global minimiser $x^*$ of $f$ is fixed throughout.

Given positive step sizes $\alpha_k>0$ and $\gamma_0>0$, the implicit scheme (72)–(73) relates consecutive iterates $(x_k,v_k,\gamma_k)$ and $(x_{k+1},v_{k+1},\gamma_{k+1})$ by

$$\frac{x_{k+1}-x_k}{\alpha_k}=v_{k+1}-x_{k+1},\qquad \frac{v_{k+1}-v_k}{\alpha_k}=\frac{\mu}{\gamma_k}(x_{k+1}-v_{k+1})-\frac1{\gamma_k}\nabla f(x_{k+1}),\qquad \frac{\gamma_{k+1}-\gamma_k}{\alpha_k}=\mu-\gamma_{k+1}.$$

All right-hand sides are evaluated at the new iterate, which is what makes the scheme implicit. The Lyapunov function (74) is

$$\mathcal L_k:=f(x_k)-f(x^*)+\frac{\gamma_k}{2}\|v_k-x^*\|^2 .$$

Formalization targets

Goal: Theorem 4.1 (p. 17)

For $f\in\mathcal S^1_\mu$ with $\mu\ge0$ and every run of (72)–(73) with $\gamma_0>0$ and $\alpha_k>0$,

$$\mathcal L_{k+1}\ \le\ \frac{\mathcal L_k}{1+\alpha_k},\qquad k\in\mathbb N .$$

No condition links $\alpha_k$ to $\gamma_k$ or to any smoothness constant, and $f$ need not have a Lipschitz gradient.

Milestones (proof of Theorem 4.1, p. 17)

1. The difference identity: under (73),

$$\mathcal L_{k+1}-\mathcal L_k=f(x_{k+1})-f(x_k)+\tfrac{\alpha_k}2(\mu-\gamma_{k+1})\|v_{k+1}-x^*\|^2+\gamma_k\big(v_{k+1}-v_k,\tfrac{v_{k+1}+v_k}2-x^*\big).$$

2. The last-term identity $\gamma_k\big(v_{k+1}-v_k,\tfrac{v_{k+1}+v_k}2-x^*\big)=\gamma_k(v_{k+1}-v_k,v_{k+1}-x^*)-\tfrac{\gamma_k}2\|v_{k+1}-v_k\|^2$.
3. The substitution of the second line of (72) into $\gamma_k(v_{k+1}-v_k,v_{k+1}-x^*)$.
4. The split of the gradient term by the first line of (72).
5. The one-step inequality (75): $\mathcal L_{k+1}-\mathcal L_k\le-\alpha_k\mathcal L_{k+1}$.

A companion item states (76): each step of (72) is a proximal point step $x_{k+1}=\mathbf{prox}_{\eta_kf}(y_k)$ with explicit $\eta_k$, $y_k$ and $\gamma_{k+1}=(\gamma_k+\mu\alpha_k)/(1+\alpha_k)$.

Significance

Theorem 4.1 gives a linear rate $\mathcal L_k\le\mathcal L_0\prod_{i<k}(1+\alpha_i)^{-1}$ for every convex $f$, not only strongly convex ones, as long as $\alpha_k\ge\alpha>0$; larger steps give faster decay. It is the exact discrete counterpart of the exponential decay $\mathcal L(t)\le e^{-t}\mathcal L(0)$ of the NAG flow (Lemma 3.2 of the paper), and the semi-implicit and explicit schemes of §§5–7 are analysed by comparison with it: their proofs follow the same difference computation and must control the extra terms that the implicit scheme does not produce. Through (76), the implicit scheme is an extrapolated proximal point method, close to Güler's accelerated proximal point algorithm when $\mu=0$.

The result is proved on paper. To the knowledge of this mission it has no machine-checked proof. Formalizing it fixes the exact hypotheses of the discrete Lyapunov argument (in particular that $\gamma_k>0$ and $\mu\ge0$ are what make the dropped squares nonpositive) and provides the algebraic identities reused by the later schemes of the paper.

Difficulty

The argument is algebraic, but the bookkeeping is where it can fail. The difference $\mathcal L_{k+1}-\mathcal L_k$ produces terms in $\|v_{k+1}-x^*\|^2$, $\|x_{k+1}-x^*\|^2$, $\|x_{k+1}-v_{k+1}\|^2$, $\|v_{k+1}-v_k\|^2$ and $\|x_{k+1}-x_k\|^2$, which must cancel or carry the right sign exactly. The naive attempt, to bound $\mathcal L_{k+1}-\mathcal L_k$ by $-\alpha_k\mathcal L_k$ as in an explicit Euler analysis, fails: the implicit scheme controls the decrease by the new value $\mathcal L_{k+1}$, and the contraction factor $1/(1+\alpha_k)$ comes from solving for it. The $\gamma$-terms close only because (73) is also discretized implicitly, so that $\gamma_{k+1}$, not $\gamma_k$, multiplies $\|v_{k+1}-x^*\|^2$.

Formalization scope

$V$ is a type with NormedAddCommGroup, InnerProductSpace ℝ V and CompleteSpace V. The gradient is an explicit map gradf : V → V with HasGradientAt f (gradf x) x at every point, and continuity of gradf encodes "continuously differentiable". The class $\mathcal S^1_\mu$, the Lyapunov function, one step of (72)–(73), and a run are definitions; the scheme's equations are kept in the page's difference-quotient form, as hypotheses on given sequences. Nothing asserts that the implicit system has a solution, and iterates are not defined by choice. Sequences are indexed from $0$.

Conventions committed to: $x^*$ is a global minimiser (the paper's standing assumption that $\operatorname{argmin}f$ is nonempty); $\gamma_0>0$ and $\alpha_k>0$ as on p. 16; the one-step milestones that divide by $\gamma_k$ assume $\gamma_k>0$, which along a run follows from (73). The prox in (76) is stated as membership in the argmin through the published predicate GoldenRatioVI.Shared.IsProxPoint.

A trivializing formalization is ruled out: the goal is stated for sequences satisfying the scheme (72)–(73) itself, never for sequences assumed to satisfy (75), and no Lipschitz-gradient class or step-size rule from later sections is imported.

The definitions of $\mathcal S^1_\mu$ and $\mathcal L_k$ are shared with the other missions of this series and will be consolidated. Contributions welcome: proofs of the identities, of (75), of the goal, and of (76).

Selected references

• H. Luo, L. Chen, From differential equation solvers to accelerated first-order methods for convex optimization, Math. Program., 2022; arXiv:1909.03145v4. https://arxiv.org/abs/1909.03145
• W. Su, S. Boyd, E. J. Candès, A differential equation for modeling Nesterov's accelerated gradient method: theory and insights, J. Mach. Learn. Res. 17, 2016. https://arxiv.org/abs/1503.01243
• O. Güler, New proximal point algorithms for convex minimization, SIAM J. Optim. 2(4), 649–664, 1992. https://doi.org/10.1137/0802032
• R. T. Rockafellar, Monotone operators and the proximal point algorithm, SIAM J. Control Optim. 14(5), 877–898, 1976. https://doi.org/10.1137/0314056
• Y. Nesterov, A method for solving the convex programming problem with convergence rate O(1/k²), Soviet Math. Dokl. 27, 372–376, 1983.

Motivation

Douglas–Rachford splitting (DRS) minimizes a sum $\varphi_1+\varphi_2$ by alternating between the proximal mappings of the two terms. It is one of the most used first-order methods in signal processing, statistics and control, and, through its dual interpretation, it is the engine behind the alternating direction method of multipliers (ADMM). For convex problems its convergence theory is classical. In nonconvex problems — sparse regression with $\ell_0$ penalties, rank constraints, phase retrieval, feasibility problems with nonconvex sets — DRS is used routinely, but its theory is recent and depends on how the stepsize $\gamma$ and the relaxation $\lambda$ are chosen.

Themelis and Patrinos (SIAM J. Optim. 30(1), 2020; preprint arXiv:1709.05747v4) prove that when $\varphi_1$ is $L$-smooth and $\sigma$-hypoconvex and $\varphi_2$ is merely proper and lower semicontinuous, DRS converges subsequentially to stationary points for stepsizes $\gamma<1/L$ and relaxations in an explicit range. Their §4.2 then shows that these ranges cannot be enlarged: it constructs problems on which DRS fails as soon as the stepsize or the relaxation leaves the range. This mission formalizes that tightness half.

Timeline. Lions and Mercier (1979) proved convergence of DRS for maximal monotone operators. Li and Pong (Math. Program. 2016) gave the first nonconvex analysis, for $\lambda=1$ and a stepsize bound well below $1/L$. Themelis and Patrinos (2017–2020) introduced the Douglas–Rachford envelope, extended the analysis to all relaxations, reached the bound $\gamma<1/L$, and proved with the examples of §4.2 that the bounds are tight.

Setting

Write $\overline{\mathbb R}=\mathbb R\cup\{+\infty\}$. A function $h:\mathbb R^p\to\overline{\mathbb R}$ is proper if it is never $-\infty$ and finite somewhere. A differentiable $h:\mathbb R^p\to\mathbb R$ is $L$-smooth if $\nabla h$ is $L$-Lipschitz, and $\sigma$-hypoconvex if $h-\tfrac\sigma2\|\cdot\|^2$ is convex; an $L$-smooth function is $\sigma$-hypoconvex for some $\sigma\in[-L,L]$, and $\sigma>0$ is strong convexity.

The proximal mapping of $h$ with parameter $\gamma>0$ is the set

$$\operatorname{prox}_{\gamma h}(x)=\operatorname*{arg\,min}_{w}\Big\{h(w)+\tfrac1{2\gamma}\|w-x\|^2\Big\},$$

which may be empty or contain several points. Given $\varphi_1,\varphi_2$, one DR-iteration with stepsize $\gamma>0$ and relaxation $\lambda>0$ maps $s$ to $s^+$ through

$$u\in\operatorname{prox}_{\gamma\varphi_1}(s),\qquad v\in\operatorname{prox}_{\gamma\varphi_2}(2u-s),\qquad s^+=s+\lambda(v-u).$$

A DRS run is a sequence $(s^k,u^k,v^k)_{k\in\mathbb N}$ obtained by iterating this step with any selection from the proximal sets. Its residual is $\|s^k-s^{k+1}\|=\lambda\|u^k-v^k\|$. A point $s$ is a fixed point if some admissible step gives $s^+=s$.

The paper's standing assumption (Assumption I) is: $\varphi_1$ is $L$-smooth and $\sigma$-hypoconvex, $\varphi_2$ is proper and lower semicontinuous, and $\operatorname{arg\,min}(\varphi_1+\varphi_2)\ne\emptyset$. The positive result needs $\gamma<1/L$; Theorem 4.7 asks what happens beyond it.

Formalization targets

Goal: Theorem 4.7, necessity of $\gamma<1/L$

For every $p\ge1$, $L>0$ and $\sigma\in[-L,L]$ there are $\varphi_1:\mathbb R^p\to\mathbb R$ and $\varphi_2:\mathbb R^p\to\overline{\mathbb R}$ satisfying Assumption I such that

$$\forall\,\gamma>\tfrac1L,\ \forall\,\lambda>0,\ \forall\text{ DRS runs }(s^k,u^k,v^k):\qquad \|s^k-s^{k+1}\|\not\to0 .$$

The pair may depend on $p$, $L$ and $\sigma$ only, never on $\gamma$, $\lambda$ or the starting point.

Milestones on the paper's one-dimensional example

1. The function (4.12), with $\varphi_2=\delta_{\{\pm1\}}$, satisfies Assumption I, and $\pm1$ are global minimizers.
2. (4.13): $\operatorname{prox}_{\gamma\varphi_1}$ is single-valued everywhere iff $\gamma<1/[\sigma]_-$, with an explicit two-branch formula, and $\operatorname{prox}_{\gamma\delta_{\{\pm1\}}}=\operatorname{sgn}$ with $\operatorname{sgn}(0)=\{\pm1\}$.
3. The case analysis of $u^k-v^k$ along any run with $1/L<\gamma<1/[\sigma]_-$.

Companion: Theorem 4.8, necessity of $\lambda<2(1+\gamma\sigma)$

For every $p$, $L>0$ and $\sigma\in[-L,L]$ there are $\varphi_1$, $\varphi_2$ with $\varphi_2$ additionally strongly convex such that, for every $0<\gamma<1/L$, every $\lambda>2(1+\gamma\sigma)$ and every run, either $s^0$ is a fixed point or the residual does not vanish. A further milestone is the expansion identity $|u^{k+1}-p|=|1-\tfrac{\lambda}{1+\gamma\sigma}|\,|u^k-p|$ from its proof.

Significance

The result. Theorems 4.7 and 4.8, together with the paper's convergence theorem, give a sharp map of the parameters for which nonconvex DRS is guaranteed to work under Assumption I. For practitioners they say that the stepsize $\gamma<1/L$ cannot be traded for a larger one without structure on $\varphi_2$, and that over-relaxation beyond $2(1+\gamma\sigma)$ fails even when $\varphi_2$ is strongly convex. The paper's Theorem 4.9 summarizes these facts as tightness of the bound $\gamma<\min\{1/L,(2-\lambda)/(2[\sigma]_-)\}$.

Formalizing it. The results are proved on paper but not machine-checked anywhere. The formalization contributes a precise, selection-robust statement — the claim covers every run, including runs through points where a proximal set has two elements — and records two places where the printed statement and the proof do not match (see Formalization scope). Explicit proximal computations for piecewise-quadratic functions and indicators of finite sets, which the milestones need, are reusable for other counterexamples in splitting methods.

Difficulty

The claim is about every run, not one. A proximal set may contain two points, and a counterexample that works for one selection can fail for another. This is exactly what happens at $\gamma=1/L$ in the paper's example: one admissible selection stays at a minimizer for ever.

The second difficulty is the dimension: the paper's functions live on $\mathbb R$ while the theorems are stated on $\mathbb R^p$, and the construction has to be lifted so that the extra coordinates neither create fixed points nor make the residual vanish. Finally, for $\sigma<0$ and large $\gamma$ the proximal mapping of $\varphi_1$ is empty or multivalued, so the analysis of a run must not assume a closed-form update.

Formalization scope

$\mathbb R^p$ is EuclideanSpace ℝ (Fin p); $\varphi_1$ is real-valued and $\varphi_2$ is EReal-valued and never $-\infty$. The proximal mapping is a set and a run is any selection; no proximal function is chosen. Every reciprocal bound is multiplied out: $\gamma>1/L$ is 1 < γ * L, $\gamma<1/L$ is γ * L < 1, $\gamma<1/[\sigma]_-$ is γ * negPartR σ < 1, so that $[\sigma]_-=0$ imposes no constraint. The milestones are stated on $\mathbb R$, as in the paper.

Departures from the page, each disclosed in the items:

1. Theorem 4.7 is stated for $\gamma>1/L$; the page has $\gamma\ge1/L$. The printed version is false at $\gamma=1/L$ for every pair satisfying Assumption I: if $x^\star$ is a minimizer and $s^\star=x^\star+\gamma\nabla\varphi_1(x^\star)$, then $s^k=s^\star$, $u^k=v^k=x^\star$ is a DRS run with zero residual. In the paper's example ($L=1$, $\sigma=0$, $t=2$, $\gamma=1$) it is $s^k=2$, $u^k=v^k=1$.
2. $p\ge1$ is added to Theorem 4.7: on $\mathbb R^0$ every run is constant. Theorem 4.8 is stated for every $p$.
3. Theorem 4.8 is stated as printed ("unless $s^0$ is a fixed point"). The paper's own construction does not establish it: with $L=1$, $\sigma=-1$, $\gamma=\tfrac12$, $\lambda=\tfrac65$, $p=2$, $s^0=-2$, the first step reaches a fixed point although $s^0$ is not one. The theorem is existential, so the statement stands; solvers may use any pair satisfying p1–p3.
4. "Well defined" in (4.13) is read as "nonempty and single-valued at every point".
5. For $\gamma\ge1/[\sigma]_-$ the proximal sets of $\varphi_1$ may be empty, so no run exists and the claim holds with nothing to check; this matches the page, whose proof treats $1/L\le\gamma<1/[\sigma]_-$. In particular, at $\sigma=-L$ the goal is satisfied by pairs admitting no runs at all — as on the page.

A formalization in which $\varphi_2$ is improper, or in which the argmin is empty, would make the claim degenerate; properness and the argmin clause are part of the existential and rule that out. Weakening p4 to "some run" or "some starting point" is not this theorem.

Contributions welcome: proofs of the milestones (piecewise-quadratic proximal mappings, the prox of $\delta_{\{\pm1\}}$), the lift from $\mathbb R$ to $\mathbb R^p$, and the two theorems.

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. https://doi.org/10.1137/18M1163993 — preprint arXiv:1709.05747v4, https://arxiv.org/abs/1709.05747
• P.-L. Lions, B. Mercier, Splitting algorithms for the sum of two nonlinear operators, SIAM J. Numer. Anal. 16(6):964–979, 1979. https://doi.org/10.1137/0716071
• G. Li, T. K. Pong, Douglas–Rachford splitting for nonconvex optimization with application to nonconvex feasibility problems, Math. Program. 159:371–401, 2016. https://doi.org/10.1007/s10107-015-0963-5
• R. T. Rockafellar, R. J.-B. Wets, Variational Analysis, Springer, 1998. https://doi.org/10.1007/978-3-642-02431-3