Numerical Analysis

Motivation

The gradient method minimizes a differentiable function $f$ on $\mathbb{R}^n$ by moving from the current point $x_k$ along the negative gradient, $x_{k+1} = x_k - \lambda_k \nabla f(x_k)$. Every practical implementation has to choose the step $\lambda_k$. A fixed step requires knowing a Lipschitz constant of $\nabla f$; an exact line search is expensive. Larry Armijo's four-page note (Pacific J. Math. 16 (1966) 1–3) proved that any step achieving a sufficient decrease proportional to $|\nabla f(x_k)|^2$ yields convergence, and gave the halving rule that is now called the Armijo rule or backtracking line search. The rule is the default step-size strategy of gradient, quasi-Newton and Newton-type methods in textbooks such as Nocedal–Wright and Bertsekas, and its sufficient-decrease test reappears in essentially every later line-search analysis.

Timeline:

• 1944, H. B. Curry: convergence of steepest descent when each step minimizes $f$ along the ray, which needs a one-dimensional minimization per step (Quart. Appl. Math. 2).
• 1962, A. A. Goldstein: convergence of the gradient method with a rate estimate, assuming $f \in C^2$ on a bounded level set $S(x_0)$ and a known bound on the Hessian norm there (Numer. Math. 4).
• 1966, L. Armijo (received January 1964): the convergence theorem below for any step with sufficient decrease $\delta|\nabla f(x_k)|^2$, $0 < \delta \le 1/(4K)$, under a Lipschitz condition on $\nabla f$ restricted to the level set, with the fixed-step and halving-step algorithms as corollaries. The paper's §3 places these hypotheses between Curry's (weaker) and Goldstein's (stronger).

Setting

Let $E^n$ be real Euclidean $n$-space with norm $|\cdot|$, and let $f : E^n \to \mathbb{R}$ be continuous and bounded below. For a fixed $x_0 \in E^n$ the level set is $S(x_0) = \{x : f(x) \le f(x_0)\}$. Write $\nabla f(x)$ for the gradient of $f$ at $x$. "$f \in C^1$ on $S(x_0)$" means that $f$ is differentiable at every point of $S(x_0)$ and $\nabla f$ is continuous there.

• Condition III at $x_0$: $f \in C^1$ on $S(x_0)$ and there is a Lipschitz constant $K > 0$ with $|\nabla f(y) - \nabla f(x)| \le K|y - x|$ for all $x, y \in S(x_0)$.
• Condition IV at $x_0$: $f \in C^1$ on $S(x_0)$, $x^*$ is a point with $f(x^*) = \inf_{E^n} f$, and for every $r > 0$

$$m(r) = \inf\{|\nabla f(x)| : x \in S(x_0),\ |x - x^*| \ge r\} > 0,$$

with $m(r) = \infty$ when the set is empty. Condition IV forces $x^*$ to be the unique minimizer and the only stationary point in $S(x_0)$.

For $\delta > 0$ and $x \in S(x_0)$, the sufficient-decrease set of display (1) is

$$S^*(x,\delta) = \{x_\lambda : x_\lambda = x - \lambda\nabla f(x),\ \lambda > 0,\ f(x_\lambda) - f(x) \le -\delta|\nabla f(x)|^2\}.$$

The steepest descent algorithm uses $x_{k+1} = x_k - \frac{1}{2K}\nabla f(x_k)$. The modified steepest descent algorithm fixes $\alpha > 0$, sets $\alpha_m = \alpha/2^{m-1}$, and takes $x_{k+1} = x_k - \alpha_{m_k}\nabla f(x_k)$ with $m_k$ the smallest positive integer such that

$$f(x_k - \alpha_{m_k}\nabla f(x_k)) - f(x_k) \le -\tfrac12\alpha_{m_k}|\nabla f(x_k)|^2. \tag{2}$$

Formalization targets

Goal: the convergence theorem (§2, pp. 1–2)

Assume $f$ is continuous and bounded below and Conditions III and IV hold at $x_0$. If $0 < \delta \le 1/(4K)$, then for every $x \in S(x_0)$

$$\emptyset \ne S^*(x,\delta) \subseteq S(x_0),$$

and every sequence with first term $x_0$ and $x_{k+1} \in S^*(x_k,\delta)$ for all $k$ satisfies

$$x_k \to x^* \qquad (k \to \infty).$$

Milestones (the claims of the proof, p. 2)

1. Along such a sequence, $|\nabla f(x_k)| \to 0$.
2. Under Condition IV, a sequence in $S(x_0)$ with $|\nabla f| \to 0$ converges to $x^*$.

Further claims of the proof (p. 2)

• For $x \in S(x_0)$ and $0 \le \lambda \le 1/K$: $f(x_\lambda) - f(x) \le -(\lambda - \lambda^2K)|\nabla f(x)|^2$.
• $x_\lambda \in S^*(x,\delta)$ for $\lambda_1 \le \lambda \le \lambda_2$, $\lambda_i = \frac{1}{2K}[1 + (-1)^i\sqrt{1 - 4\delta K}]$.
• Along such a sequence, $x_k \in S(x_0)$ and $f(x_k)$ is nonincreasing.

Companion results

• Corollary 1: the steepest descent iterates converge to $x^*$.
• Corollary 2: the modified steepest descent iterates converge to $x^*$, for every $\alpha > 0$.
• The index $m_k$ of Corollary 2 exists; $m_k = 1$ if $\alpha \le 1/(2K)$ and $\alpha_{m_k} > 1/(4K)$ otherwise; the half-step estimate $f(x_\lambda) - f(x) \le -\frac12\lambda|\nabla f(x)|^2$ for $0 \le \lambda \le 1/(2K)$.
• The §1 remark: Condition IV implies Conditions I and II, and is equivalent to them when $S(x_0)$ is bounded.

Significance

The theorem separates the analysis of a gradient method from the choice of its step: whatever rule produces steps in $S^*(x_k, \delta)$ converges, and the two corollaries show that both a fixed step $1/(2K)$ and a halving search started from an arbitrary $\alpha$ qualify. The halving rule needs no knowledge of $K$, which is why it became the standard globalization device for descent methods. The hypotheses are local to the level set: $\nabla f$ need only be Lipschitz on $S(x_0)$, not on all of $E^n$, so functions whose gradient is only locally Lipschitz, such as polynomials of degree above two with bounded level sets, satisfy Condition III.

The result is classical and proved on the page; it has not, to our knowledge, been machine-checked in this form. Formalizing it provides a verified sufficient-decrease convergence theorem with level-set-local hypotheses, a verified Armijo backtracking rule with explicit constants $\tfrac12$, $1/(4K)$ and $1/(8K)$, and reusable definitions (sufficient-decrease set, Armijo test with halving steps) on which later line-search results can build.

Difficulty

The mean value step in the proof of the descent inequality silently uses the Lipschitz bound along the whole segment from $x$ to $x - \lambda\nabla f(x)$, but Condition III gives it only on $S(x_0)$. One has to show that the segment cannot leave the closed set $S(x_0)$ for $0 \le \lambda \le 1/K$; replacing the hypothesis by a global Lipschitz condition is not acceptable, because it changes the theorem. The final step needs Condition IV in its exact form: gradients tending to zero along the iterates do not by themselves give convergence of the iterates, and the positive lower bound $m(r)$ away from $x^*$ is what turns stationarity into convergence. Corollary 2 additionally requires a termination argument for the halving search and a uniform lower bound on accepted steps.

Formalization scope

$E^n$ is EuclideanSpace ℝ (Fin n), $|\cdot|$ is the norm and $\nabla f$ is Mathlib's gradient f. All declarations live in the namespace ArmijoGrad.Conv, with one definition file ArmijoGrad.Conv.Setting.

Committed conventions:

• The standing assumptions of §2 are binders of the goal and the corollaries: Continuous f, BddBelow (Set.range f), ConditionIII f x0 K, ConditionIV f x0 xstar. Milestones keep only the assumptions their claim uses, which makes them more general; nothing is added.
• "$f \in C^1$ on $S(x_0)$" is differentiability at every point of $S(x_0)$ plus continuity of gradient f on $S(x_0)$. The Lipschitz bound of Condition III holds on $S(x_0)$ only; $K > 0$ is part of Condition III.
• Condition IV carries its minimizer $x^*$ and states $m(r) > 0$ as the existence of a positive lower bound for $|\nabla f|$ on $\{x \in S(x_0) : |x - x^*| \ge r\}$; this matches the convention $m(r) = \infty$ on the empty set. No real infimum is used.
• Sequences are ℕ → EuclideanSpace ℝ (Fin n) with first term equal to $x_0$, the base point of $S(x_0)$, as in the paper's notation.
• $\lambda > 0$ in (1) is strict. $\alpha_m = \alpha/2^{m-1}$ is used only with $m \ge 1$, and $m_k$ is required to be the smallest such index.

A trivializing formalization is ruled out: the goal assumes nothing about $\lambda_1, \lambda_2$, the descent inequality or $|\nabla f(x_k)| \to 0$, Condition IV cannot hold vacuously because it names a minimizer, and all hypotheses are satisfiable by nontrivial data: $f(x) = |x|^2/2$ satisfies Conditions III and IV with $K = 1$ and $x^* = 0$, and both algorithms have runs for it.

A complete development needs the one-dimensional mean value inequality along a segment, a continuity argument keeping the segment in $S(x_0)$, the telescoping bound for monotone sequences bounded below, and the halving-search termination. The definitions file and the descent and Armijo-index lemmas are reusable for later line-search missions. Contributions are welcome on every milestone; the descent inequality is the natural first target.

Selected references

• L. Armijo, Minimization of functions having Lipschitz continuous first partial derivatives, Pacific Journal of Mathematics 16(1), 1–3, 1966. https://doi.org/10.2140/pjm.1966.16.1
• H. B. Curry, The method of steepest descent for non-linear minimization problems, Quarterly of Applied Mathematics 2, 1944. https://doi.org/10.1090/qam/10667
• A. A. Goldstein, Cauchy's method of minimization, Numerische Mathematik 4, 146–150, 1962. https://doi.org/10.1007/BF01386306
• J. Nocedal and S. J. Wright, Numerical Optimization, 2nd ed., Springer, 2006. https://doi.org/10.1007/978-0-387-40065-5

Motivation

Gradient methods are the workhorse of large-scale nonlinear optimization and of the training of statistical models. In practice the direction actually used is rarely the exact negative gradient: it may be scaled, computed incrementally one data component at a time, or perturbed by approximation error. The classical convergence theory of such methods often assumes that the iterates stay bounded, that the objective is bounded below, or that the errors vanish at a prescribed rate, and these assumptions must then be checked separately for each method.

Bertsekas and Tsitsiklis (2000) proved convergence results for gradient methods with errors that need none of these assumptions. Their deterministic result (Proposition 1) allows a general descent direction together with an error whose size is proportional to the stepsize, and concludes that either the objective values diverge to $-\infty$ or they converge and the gradients tend to zero. It applies, among others, to the incremental gradient method for a sum of functions (Proposition 2 of the same paper), which underlies backpropagation-style training. The stochastic counterpart (Proposition 3, zero-mean errors) is the subject of a companion mission.

Setting

Throughout, $f:\mathbb R^n\to\mathbb R$ is a continuously differentiable function whose gradient is Lipschitz continuous: for some constant $L$,

$$\|\nabla f(x)-\nabla f(\bar x)\|\le L\|x-\bar x\|\qquad\forall x,\bar x\in\mathbb R^n. \tag{2.1}$$

Here $\|\cdot\|$ is the Euclidean norm and $x'y$ the standard inner product. The gradient method with errors generates a sequence of iterates

$$x_{t+1}=x_t+\gamma_t(s_t+w_t),\qquad t=0,1,\dots,$$

where $\gamma_t>0$ is the stepsize, $s_t$ is a descent direction and $w_t$ is an error vector. Nothing is assumed about how $s_t$ and $w_t$ are produced, beyond the two conditions below, which hold for some positive scalars $c_1,c_2,p,q$ and every $t$:

$$c_1\|\nabla f(x_t)\|^2\le-\nabla f(x_t)'s_t,\qquad\|s_t\|\le c_2\bigl(1+\|\nabla f(x_t)\|\bigr), \tag{2.2}$$ $$\|w_t\|\le\gamma_t\bigl(q+p\|\nabla f(x_t)\|\bigr). \tag{2.3}$$

The stepsizes are diminishing in the standard sense:

$$\sum_{t=0}^\infty\gamma_t=\infty,\qquad\sum_{t=0}^\infty\gamma_t^2<\infty.$$

A stationary point of $f$ is a point $\bar x$ with $\nabla f(\bar x)=0$; a limit point of $(x_t)$ is the limit of some subsequence.

Formalization targets

Goal: Proposition 1 (p. 630)

Under (2.1), (2.2), (2.3) and the stepsize conditions, either

$$f(x_t)\to-\infty,$$

or else $f(x_t)$ converges to a finite value and

$$\lim_{t\to\infty}\nabla f(x_t)=0.$$

Furthermore, every limit point of $(x_t)$ is a stationary point of $f$.

Milestones

The milestones follow the paper's own argument, in order.

1. Lemma 1 (p. 629). For real sequences with $W_t\ge0$, $Y_{t+1}\le Y_t-W_t+Z_t$ and $\sum_{t=0}^T Z_t$ convergent, either $Y_t\to-\infty$, or $Y_t$ converges and $\sum_t W_t<\infty$.
2. (2.4) (p. 630). Under (2.1), $f(x+z)\le f(x)+z'\nabla f(x)+\tfrac L2\|z\|^2$ for all $x,z$.
3. (2.5) (p. 631). For some $\beta_1,\beta_2>0$ and all sufficiently large $t$, $f(x_{t+1})\le f(x_t)-\gamma_t\beta_1\|\nabla f(x_t)\|^2+\gamma_t^2\beta_2$.
4. (2.6) (p. 631). Either $f(x_t)\to-\infty$, or $f(x_t)$ converges and $\sum_t\gamma_t\|\nabla f(x_t)\|^2<\infty$.
5. After (2.6) (p. 631). If $f(x_t)\not\to-\infty$, then $\liminf_{t\to\infty}\|\nabla f(x_t)\|=0$.

Significance

The result. Proposition 1 separates two concerns that are usually entangled: what the method guarantees, and what must be known about $f$. It concludes stationarity of all limit points and convergence of the gradients to zero without assuming that the iterates are bounded or that $f$ is bounded below; when $f$ is bounded below the first alternative is excluded and $\nabla f(x_t)\to0$ follows outright. Because $s_t$ need not be the negative gradient and $w_t$ need not vanish faster than the stepsize, the result covers scaled gradient methods, incremental gradient methods for sums of functions, and gradient methods with deterministic approximation error. The descent inequality (2.4) and the deterministic supermartingale-type Lemma 1 are standard tools that recur throughout optimization theory.

Formalizing it. The proposition is proved in the paper; it has not been machine-checked. A formal proof would give a reusable, verified convergence theorem for a broad class of first-order methods on $\mathbb R^n$, together with a formal descent lemma for functions with Lipschitz gradient, which Mathlib does not currently state in this form, and a deterministic Robbins–Siegmund-type lemma for sequences.

Difficulty

The summability estimate (2.6) gives only $\liminf\|\nabla f(x_t)\|=0$. Passing to $\lim\nabla f(x_t)=0$ is the main step: the obvious argument (a summable series $\sum_t\gamma_t\|\nabla f(x_t)\|^2$ with $\sum_t\gamma_t=\infty$ forces the gradient norms to zero) is false in general, because a nonnegative sequence with these two properties may still have infinitely many large terms. What is missing is a bound on how far $\|\nabla f(x_t)\|$ can travel while the stepsizes are small, and only (2.1) and (2.2)–(2.3) together supply it. A second difficulty is that there is no boundedness of the iterates: every estimate must hold globally, and the error $w_t$ is controlled only relative to $\|\nabla f(x_t)\|$, which may be unbounded along the sequence.

Formalization scope

• $\mathbb R^n$ is EuclideanSpace ℝ (Fin n), $x'y$ is the real inner product ⟪x, y⟫_ℝ, and $\nabla f$ is Mathlib's gradient f; the hypothesis ContDiff ℝ 1 f makes it the true gradient. The paper's statement is for $\mathbb R^n$ and the formalization does not generalize to Hilbert spaces.
• The standing assumption (2.1) of §2 is part of every statement about $f$, as LipschitzWith L (gradient f) with L : ℝ≥0; this is equivalent to (2.1) for some real constant.
• The sequences $x_t,s_t,w_t$ and $\gamma_t$ are arbitrary data indexed from $t=0$, constrained only by the recursion and by (2.2), (2.3), $\gamma_t>0$ (all four constants $c_1,c_2,p,q$ are positive, as printed).
• $\sum_t\gamma_t=\infty$ is divergence of the partial sums to $+\infty$; $\sum_t\gamma_t^2<\infty$ is summability of nonnegative terms. In Lemma 1 the convergence of $\sum_t Z_t$ is convergence of the partial sums, not absolute convergence, since $Z_t$ may change sign.
• $\lim\inf$ is written out as "for every $\epsilon>0$, infinitely often below $\epsilon$", avoiding junk values of a lim inf of an unbounded sequence. Limit points are cluster points of the sequence.
• "Every limit point is stationary" is a separate conjunct, outside the dichotomy, exactly as on the page.
• The constants $\beta_1,\beta_2$ in (2.5) are existential. The milestones (2.5) and (2.6) retain the stepsize hypotheses of Proposition 1, where the paper derives them.
• A trivializing formalization would drop the "$-\infty$" alternative or require $f$ bounded below; neither is done. Instances satisfying all hypotheses with $f(x_t)\to-\infty$ and $\nabla f\not\to0$ exist (a linear $f$), so the first alternative is genuinely needed.
• Out of scope: Proposition 2 (the incremental gradient method of §3), which is a corollary of the goal, and the stochastic results of §4–§5.

Contributions welcome: proofs of the descent lemma and Lemma 1 (both reusable well beyond this mission), and of the steps (2.5)–(2.6) and the excursion argument.

Selected references

• D. P. Bertsekas and J. N. Tsitsiklis, Gradient Convergence in Gradient Methods with Errors, SIAM Journal on Optimization 10(3):627–642, 2000. https://doi.org/10.1137/S1052623497331063
• D. P. Bertsekas, Nonlinear Programming, 2nd ed., Athena Scientific, 1999.
• H. Robbins and D. Siegmund, A convergence theorem for non negative almost supermartingales and some applications, in Optimizing Methods in Statistics, Academic Press, 1971, 233–257. https://doi.org/10.1016/B978-0-12-604550-5.50015-8

Motivation

Many problems in signal processing, statistics and machine learning take the form

$$\min_{x\in\mathcal H}\ \Psi(x)+\Phi(x),$$

where $\Phi$ is smooth and convex (a data-fit term) and $\Psi$ is convex but possibly nonsmooth or infinite-valued (an $\ell^1$ penalty, the indicator function of a constraint set). The forward-backward method alternates a gradient step on $\Phi$ with a proximal step on $\Psi$ and reduces the objective gap at rate $O(k^{-1})$ after $k$ iterations. Combining it with Nesterov's extrapolation scheme gives the accelerated forward-backward method, best known as FISTA, which improves the guaranteed rate to $O(k^{-2})$. FISTA and its variants are standard solvers for sparse recovery and image reconstruction.

Timeline:

• 1983. Nesterov introduces the extrapolation scheme for smooth convex minimization, with an $O(k^{-2})$ rate for function values (Nesterov 1983).
• 2009. Beck and Teboulle extend it to the composite problem above (FISTA), with the $O(k^{-2})$ rate (doi:10.1137/080716542).
• 2014. Su, Boyd and Candès read the scheme as a discretization of the ODE $\ddot x+\frac{\alpha}{t}\dot x+\nabla\Theta(x)=0$ (arXiv:1503.01243).
• 2014–2015. Chambolle and Dossal (doi:10.1007/s10957-015-0746-4) and, independently, Attouch, Chbani, Peypouquet and Redont (arXiv:1507.04782) prove weak convergence of the iterates for the variant with parameter $\alpha>3$. Before that, convergence of the iterates had been open for about two decades.
• 2015. May shows that for $\alpha>3$ the continuous-time gap is $o(t^{-2})$ (arXiv:1509.05598).
• 2016. Attouch and Peypouquet prove the discrete analogue: for $\alpha>3$ the function gap of the algorithm is $o(k^{-2})$ and the velocity is $o(k^{-1})$ (arXiv:1510.08740, SIAM J. Optim. 26(3):1824–1834).

Setting

Let $\mathcal H$ be a real Hilbert space with inner product $\langle\cdot,\cdot\rangle$ and norm $\|\cdot\|$. Let:

1. $\Psi:\mathcal H\to\mathbb R\cup\{+\infty\}$ be proper (finite somewhere), lower semicontinuous and convex;
2. $\Phi:\mathcal H\to\mathbb R$ be convex and continuously differentiable, with gradient $\nabla\Phi$ Lipschitz continuous with constant $L$;
3. $\Theta=\Psi+\Phi$, and assume the set of minimizers $S=\operatorname{argmin}\Theta$ is nonempty.

For $s>0$, the proximal map $\operatorname{prox}_{s\Psi}(z)$ is the unique minimizer of $u\mapsto\Psi(u)+\frac1{2s}\|u-z\|^2$. Given $\alpha>0$, a step size $0<s<1/L$ and starting points $x_0,x_1$, algorithm (2) generates

$$y_k=x_k+\frac{k-1}{k+\alpha-1}(x_k-x_{k-1}),\qquad x_{k+1}=\operatorname{prox}_{s\Psi}\big(y_k-s\nabla\Phi(y_k)\big),\qquad k\ge1.$$

The common choice is $\alpha=3$; this mission concerns $\alpha>3$.

The proofs use the operator $G_s(y)=\frac1s\big(y-\operatorname{prox}_{s\Psi}(y-s\nabla\Phi(y))\big)$, the sequence $z_k=x_k+\frac{k-1}{\alpha-1}(x_k-x_{k-1})$, a minimizer $x^*$, and the quantity

$$\mathcal E(k)=\frac{2s}{\alpha-1}(k+\alpha-2)^2\big(\Theta(x_k)-\Theta(x^*)\big)+(\alpha-1)\|z_k-x^*\|^2 .$$

Write $\theta_k=\Theta(x_k)-\Theta(x^*)$ and $d_k=\frac1{2s}\|x_{k+1}-x_k\|^2$.

Formalization targets

Goal: Theorem 1 (p. 2)

For $\alpha>3$ and $0<s<1/L$,

$$\lim_{k\to\infty}k^2\big(\Theta(x_k)-\min\Theta\big)=0\qquad\text{and}\qquad\lim_{k\to\infty}k\,\|x_{k+1}-x_k\|=0.$$

The statement fixes no constants, so it is unaffected by later improvements of the explicit bounds below.

Milestones (in attack order)

1. (9), p. 2. If $sL\le1$, then for all $x,y$: $\Theta(y-sG_s(y))\le\Theta(x)+\langle G_s(y),y-x\rangle-\frac s2\|G_s(y)\|^2$.
2. (13), p. 3. For $\alpha\ge3$ and $k\ge1$: $\mathcal E(k+1)+2s\frac{\alpha-3}{\alpha-1}k\,\theta_k\le\mathcal E(k)$.
3. Fact 1, p. 3. $(\mathcal E(k))$ is nonincreasing and has a finite limit.
4. Fact 2, p. 3. $\theta_k\le\frac{(\alpha-1)\mathcal E(1)}{2s(k+\alpha-2)^2}$ and $\|z_k-x^*\|^2\le\frac{\mathcal E(1)}{\alpha-1}$ for $k\ge1$.
5. Fact 3, p. 3. For $\alpha>3$: $\sum_{k\ge1}k\,\theta_k\le\frac{(\alpha-1)\mathcal E(1)}{2s(\alpha-3)}$.
6. (14), p. 4. $\Theta(x_{k+1})+d_k\le\Theta(x_k)+\frac{(k-1)^2}{(k+\alpha-1)^2}d_{k-1}$ for $k\ge1$.
7. Fact 4, p. 4. For $\alpha>3$: $\sum_{k\ge1}k\,d_k\le\frac{\alpha(3\alpha-5)\mathcal E(1)}{4s(\alpha-1)(\alpha-3)}$.
8. Lemma 2, p. 4. For $\alpha>3$, $\lim_k\big[k^2d_k+(k+1)^2\theta_{k+1}\big]$ exists and is finite.

Significance

The result. The $O(k^{-2})$ rate of FISTA is often quoted as optimal for first-order methods. Theorem 1 shows that for $\alpha>3$ the worst-case rate along every single run is strictly better, $o(k^{-2})$, and that the steps $\|x_{k+1}-x_k\|$ decay faster than $1/k$. No better power is possible: by Attouch et al., Example 2.13 there is no $p>2$ with an $O(k^{-p})$ rate for every $\Phi$ and $\Psi$. The intermediate estimates (Facts 2–4) give explicit, quantitative bounds that are reused in the analysis of inexact and perturbed variants (Theorem 4 of the paper) and in mission II of this series, which proves weak convergence of the iterates.

Formalizing it. The result is proved on paper. As far as we know, no machine-checked proof of the $O(k^{-2})$ rate of FISTA in this Hilbert-space, extended-valued setting exists in Mathlib or on this platform, and neither does the $o(k^{-2})$ refinement. A complete development would provide reusable statements about proximal-gradient steps for functions valued in $\mathbb R\cup\{+\infty\}$. The page also has a factor slip in Fact 4 and Lemma 2 (see Formalization scope); checking it mechanically is one of the outputs.

Difficulty

The $O(k^{-2})$ bound follows from the monotonicity of $\mathcal E$ alone. That argument cannot give $o(k^{-2})$: it controls $k^2\theta_k$ only by the constant $\mathcal E(1)$, and summability of $k\theta_k$ (Fact 3) gives a decay of $k^2\theta_k$ only along a subsequence, not of the whole sequence. The missing ingredient is convergence of a weighted combination of function gaps and velocities, which is Lemma 2. The bookkeeping is delicate: every step carries explicit coefficients in $k$ and $\alpha$, and the inequalities must hold in the extended reals because $\Psi$ may be $+\infty$ at the starting point.

Formalization scope

• Space and functions. $\mathcal H$ is a real inner product space that is complete. $\Psi$ takes values in EReal and satisfies the published predicate IsProperClosedConvex (never $-\infty$, finite somewhere, lower semicontinuous, convex epigraph). $\Phi:\mathcal H\to\mathbb R$ is ConvexOn and ContDiff ℝ 1, and gradient Φ is LipschitzWith L for some L : ℝ≥0.
• Step size. $0<s<1/L$ is written 0 < s and s * L < 1, so that $L=0$ is allowed (s < 1 / L would be unsatisfiable when $L=0$). Display (9) uses the page's $s\le1/L$, written s * L ≤ 1.
• Proximal map. It is a map P with the published predicate IsProx s Ψ P, which determines $P=\operatorname{prox}_{s\Psi}$ for proper closed convex $\Psi$ and $s>0$.
• The run. It is required to follow (2) for every $k\ge1$, with $x_0,x_1$ arbitrary; the coefficient of $x_0$ vanishes at $k=1$.
• Extended reals. $\Theta$, $\mathcal E$ and the function-value statements live in EReal. No value is ever converted to $\mathbb R$ with toReal, which would turn $+\infty$ into $0$. "The limit exists" (Fact 1, Lemma 2) means convergence to a real number, because in EReal every monotone sequence converges. Infinite series of nonnegative terms are stated as bounds on every partial sum. $\min\Theta$ in the goal is ⨅ y, Θ y together with the hypothesis that a minimizer exists.
• Corrected slips. Fact 2 is stated with $\mathcal E(1)$ for $k\ge1$; the page writes $\mathcal E(0)$ for $k\ge0$, which needs an iterate $x_{-1}$. Fact 4 is stated for $\sum k\,d_k$; the page prints $\sum k\|x_{k+1}-x_k\|^2$ with the same constant, which is off by the factor $2s$ and false in general. Lemma 2 is stated for the bracket $k^2d_k+(k+1)^2\theta_{k+1}$ of its proof (16).
• Ruled out. The goal assumes nothing about $\mathcal E$, $z_k$, $\theta_k$ or $d_k$. A statement of the first limit for toReal values, or with an extended-real "limit exists", would be trivially weaker and is not what is asked.
• Welcome contributions. Proofs of any milestone; general facts about proximal maps of EReal-valued convex functions and about the descent property of $L$-smooth convex functions, which are reusable well beyond this mission.

Selected references

• H. Attouch, J. Peypouquet, The Rate of Convergence of Nesterov's Accelerated Forward-Backward Method is Actually Faster than $1/k^2$, SIAM J. Optim. 26(3):1824–1834, 2016. arXiv:1510.08740v4, doi:10.1137/15M1046095
• H. Attouch, Z. Chbani, J. Peypouquet, P. Redont, Fast convergence of inertial dynamics and algorithms with asymptotic vanishing viscosity, Math. Program. 168:123–175, 2018. arXiv:1507.04782
• A. Beck, M. Teboulle, A fast iterative shrinkage-thresholding algorithm for linear inverse problems, SIAM J. Imaging Sci. 2(1):183–202, 2009. doi:10.1137/080716542
• A. Chambolle, C. Dossal, On the convergence of the iterates of the "fast iterative shrinkage/thresholding algorithm", J. Optim. Theory Appl. 166:968–982, 2015. doi:10.1007/s10957-015-0746-4
• R. May, Asymptotic for a second order evolution equation with convex potential and vanishing damping term, 2015. arXiv:1509.05598
• Y. Nesterov, A method of solving a convex programming problem with convergence rate $O(1/k^2)$, Soviet Math. Dokl. 27:372–376, 1983. mathnet
• 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(153):1–43, 2016. arXiv:1503.01243

Motivation

Solving a square system of nonlinear equations $\psi_i(x) = 0$, $i = 1, \dots, n$, $x \in \mathbb{R}^n$, is one of the oldest tasks of numerical analysis. It appears whenever optimality conditions, equilibrium conditions or discretized differential equations have to be solved. N. Z. Shor's Minimization Methods for Non-Differentiable Functions (Springer 1985, translated by K. C. Kiwiel and A. Ruszczyński) treats such a system as the nonsmooth minimization problem

$$\min_{x \in E_n} f(x), \qquad f(x) = \max_{1 \le i \le n} |\psi_i(x)|,$$

whose optimal value is $0$ exactly when the system is consistent. Section 3.5 of the book applies the subgradient method with space dilation to this problem. At each iteration that method needs the gradient of one of the $\psi_i$, not all $n$ of them. Taking the dilation coefficient to its extreme value ($\beta = 0$) turns the method into a gradient orthogonalization method. The book shows that this method converges at a quadratic rate near a regular solution, like Newton's method, although it never forms or factorizes the Jacobian.

This mission is the eighth in a series formalizing Shor's book. It covers Sections 3.5–3.6 (printed pp. 62–71).

Setting

$E_n$ is the $n$-dimensional Euclidean space with inner product $(x, y)$ and norm $\|x\|$. The data are functions $\psi_1, \dots, \psi_n : E_n \to \mathbb{R}$ with gradients $g_{\psi_i}(x)$, and the max-residual function is $f(x) = \max_i |\psi_i(x)|$ (3.24). The Jacobian is the matrix $J(x) = \{\partial\psi_i/\partial t_j\}_{i,j=1}^n$ of partial derivatives with respect to the coordinates $x = \{t_1, \dots, t_n\}$.

An almost-gradient of a function $f$ at $x_0$ is an accumulation point of a sequence of gradients $\nabla f(x_k)$, where $x_k \to x_0$ and $f$ is differentiable at every $x_k$ (p. 18).

The regular case (p. 63) is the situation in which $f^* = \min f = 0$ is attained at a point $x^*$, the $\psi_i$ are continuously differentiable near $x^*$, and $J(x^*)$ is nonsingular.

The gradient orthogonalization method (3.27) runs in stages of $n$ steps. A stage starting at $x_0$ produces $x_1, \dots, x_n$ as follows. At step $k+1$, $0 \le k < n$, the gradient $g_{\psi_{k+1}}(x_k)$ is projected on the subspace orthogonal to the gradients $g_{\psi_1}(x_0), \dots, g_{\psi_k}(x_{k-1})$ used earlier in the stage. Call the result $\varphi_{k+1}$. The step is

$$x_{k+1} = x_k - \frac{\psi_{k+1}(x_k)}{\|\varphi_{k+1}\|^2}\, \varphi_{k+1}.$$

For $k = 0$ there is nothing to project, so $\varphi_1 = g_{\psi_1}(x_0)$. The next stage starts at $x_n$.

Formalization targets

Goal: Theorem 3.9, one stage squares the error

Let $x^* = 0$ solve the system, let the $\psi_i$ be continuously differentiable with gradients Lipschitz with constant $L$ on a ball $S_\delta = \{x : \|x\| < \delta\}$ (3.28), and let $g_{\psi_1}(0), \dots, g_{\psi_n}(0)$ be linearly independent. Then there exist $\varepsilon > 0$ and $c > 0$ such that every stage with $\|x_0\| \le \varepsilon$ satisfies

$$\|x_n\| \le c\, \|x_0\|^2 .$$

Milestones inside the proof

The proof of Theorem 3.9 passes through these displayed results, stated here in attack order:

• Eq. (3.29): near the solution, the orthogonalized gradients satisfy $\|\varphi_{k+1}\| > b > 0$, with $b$ uniform over stages.
• Eq. (3.33): after step $k+1$, $|\psi_{k+1}(x_{k+1})| \le (L/b^2)\, \psi_{k+1}^2(x_k)$.
• Lemma 3.1: a step of length $h$ in a unit direction orthogonal to $g_{\psi_i}(x_1)$ changes $\psi_i$ by at most $L h (h + \|x_1 - x_2\|)$.
• Eqs. (3.36)–(3.37): within a stage, $\max_k \|x_k\| \le c_1 \|x_0\|$, and at the end point $f(x_n) \ge c_2 \|x_n\|$.

Companion result: Theorem 3.8

In the regular case, for every $\delta > 0$ there is a neighborhood of $x^*$ in which every almost-gradient satisfies

$$(1-\delta) f(x) \le (g_f(x), x - x^*) \le (1+\delta) f(x).$$

Significance

Theorem 3.9 says that a method which evaluates one equation and one gradient per step, and never solves a linear system, has the local quadratic rate of Newton's method on a regular system. The stages combine: once one stage starts close enough to the solution, the errors of later stages satisfy $\|x_0^{(r+1)}\| \le c\|x_0^{(r)}\|^2$. Theorem 3.8 describes the geometry of $f$ near a regular solution. The ratio of $(g_f(x), x - x^*)$ to $f(x)$ tends to $1$, which is what allows the dilation parameters of the book's Theorem 3.3 to approach their extreme values near the solution. Both results connect the nonsmooth space-dilation methods of Chapter 3 to classical Newton-type theory.

All results of this mission are proved in the book. None of them has, to our knowledge, a machine-checked proof: the platform has local quadratic convergence results for cubic-regularized Newton, proximal Newton and semismooth Newton, but not for this method. The work is formalizing Shor's proof. This includes building the Gram-determinant argument behind (3.29) and the second-order estimates, which are reusable for other projection-based solvers of Kaczmarz type.

Difficulty

A naive argument analyses each step as a Newton step for one equation and stops there. That shows each equation $\psi_{k+1}$ is small right after its own step, but the theorem needs all equations to be small at the end of the stage. Later steps move the point and can undo the progress on earlier equations. The step directions are orthogonal to the earlier gradients, but those gradients were taken at earlier points, not at the current one. Controlling this drift needs Lemma 3.1 together with the bound (3.36), which says that the whole stage stays within a multiple of $\|x_0\|$.

A second difficulty is that the method is only defined near the solution. The step divides by $\|\varphi_{k+1}\|^2$, and a uniform positive lower bound on these norms comes from the nonvanishing of a Gram determinant evaluated at $n$ different points. The final passage from small residuals to a small error needs the growth bound (3.37), which uses the nonsingularity of the Jacobian once more.

Formalization scope

• Representation. $E_n$ is EuclideanSpace ℝ (Fin n). The equations are a family ψ : Fin n → E_n → ℝ indexed from $0$, so ψ k is the book's $\psi_{k+1}$. Gradients are Mathlib's gradient, and $S_\delta$ is the open ball Metric.ball 0 δ. Continuous differentiability is ContDiff ℝ 1; in Theorem 3.8 it is ContDiffOn ℝ 1 on a neighborhood of $x^*$.
• The stage. A stage is a predicate IsOrthStage ψ x on a sequence $x : \mathbb{N} \to E_n$, which fixes $x_1, \dots, x_n$ from $x_0$. The projection is starProjection onto the orthogonal complement of the span of the earlier gradients.
• Corrected misprint. The printed Step 1, (3.27a), divides by $\|g_{\psi_1}(x_0)\|$ rather than its square. The formalization uses the square, which is what (3.27b) and the proof's (3.31)–(3.33) require at $k = 0$.
• Division by zero. The book gives no rule for $\varphi_{k+1} = 0$; Lean's $t/0 = 0$ would leave the point unchanged. All statements apply the method only near the solution, where (3.29) bounds $\|\varphi_{k+1}\|$ away from zero.
• Quantifier order. In Theorem 3.9 and in (3.29), (3.36)–(3.37), the constants ($\varepsilon$, $c$, $\delta'$, $b$, $\varepsilon_0$, $c_1$, $c_2$) are existentials placed before the stage. They depend only on the data $\psi, \delta, L$. A statement with $c$ chosen after $x_0$ would be trivial, with $c = \|x_n\|/\|x_0\|^2$, and is ruled out.
• Nonsingularity. In Theorem 3.9 it is the linear independence of $g_{\psi_i}(0)$, the book's primary hypothesis. In Theorem 3.8 it is $\det J(x^*) \ne 0$, with $J$ defined by coordinate partial derivatives.
• Almost-gradients. They are quantified universally: Theorem 3.8 holds for every almost-gradient, not for one chosen one.

Not formalized: Theorem 3.10 (pp. 70–71), the $n$-step quadratic convergence of the $\beta = 0$ r-algorithm with resetting. Its directional minimization rule is not pinned down on the page. The algorithm says "determined by minimizing", while the proof uses "the smallest positive root" of a stationarity equation. A contribution that fixes this rule and states the theorem is welcome as a follow-up.

Contributions of any of the five milestones are welcome independently. The Gram-determinant lower bound (3.29) and Lemma 3.1 are general facts about orthogonalization and Lipschitz gradients, and are useful beyond this mission.

Selected references

• N. Z. Shor, Minimization Methods for Non-Differentiable Functions, Springer Series in Computational Mathematics 3, Springer, 1985, §§3.5–3.6, pp. 62–71. https://doi.org/10.1007/978-3-642-82118-9
• J. M. Ortega and W. C. Rheinboldt, Iterative Solution of Nonlinear Equations in Several Variables, Academic Press, 1970. https://doi.org/10.1137/1.9780898719468

Motivation

The ellipsoid method is the algorithm that settled the polynomial-time solvability of linear programming (Khachiyan, 1979) and that underlies the equivalence of separation and optimization in combinatorial optimization (Grötschel, Lovász and Schrijver, 1981). Its origin is in nonsmooth convex optimization. In 1976 Yudin and Nemirovskii proposed a modified method of centered sections that localizes an optimum inside a sequence of ellipsoids, and in 1977 N. Z. Shor observed independently that the same scheme is a subgradient method with space dilation along the gradient, the family of methods he had developed since 1969. Section 3.8 of Shor's monograph Minimization Methods for Non-Differentiable Functions (Springer 1985) presents the method in this second form and proves its basic localization property.

Timeline:

• 1965. A. Yu. Levin proposes the method of centered sections (cuts through the center of gravity of a polyhedron); each cut removes at least a fixed fraction of the volume, but computing centers of gravity is impractical for $n > 3$.
• 1969–1972. Shor introduces subgradient methods with space dilation along the gradient (SDG methods).
• 1976. Yudin and Nemirovskii replace the polyhedron by a minimal ellipsoid containing a half-ellipsoid, obtaining a geometric volume decrease depending only on the dimension ([Yudin–Nemirovskii 1976]).
• 1977. Shor shows that the same method is an SDG algorithm with coefficient $\beta = \sqrt{(n-1)/(n+1)}$ ([Shor 1977]).
• 1979. Khachiyan applies the method to linear inequalities with integer data, obtaining the first polynomial-time algorithm for linear programming ([Khachiyan 1979]).

Setting

Let $E_n$ be $n$-dimensional Euclidean space with inner product $(x, y)$, and $n > 1$. For a unit vector $\xi$ and a number $\alpha$, the operator of space dilation along $\xi$ with coefficient $\alpha$ is $R_\alpha(\xi) = I + (\alpha - 1)\xi\xi^T$: it multiplies the component of a vector along $\xi$ by $\alpha$ and leaves the orthogonal component unchanged.

Let $g : E_n \to E_n$ be a vector field, not necessarily continuous. The problem is to find a point $x^*$ with

$$(g(x), x - x^*) \ge 0 \quad \text{for all } x \in E_n,$$

where it is known that such an $x^*$ exists in the closed ball $S(x_0, R)$ of radius $R > 0$ about a given point $x_0$. Put $\beta = \sqrt{(n-1)/(n+1)}$ and $r = n/\sqrt{n^2-1}$. The algorithm (3.57)–(3.60) starts from $x_0$, $B_0 = I_n$, $h_0 = R/(n+1)$, and at iteration $k+1$ stops if $g(x_k) = 0$, and otherwise sets

$$\xi_k = \frac{B_k^T g(x_k)}{\|B_k^T g(x_k)\|}, \quad x_{k+1} = x_k - h_k B_k \xi_k, \quad B_{k+1} = B_k R_\beta(\xi_k), \quad h_{k+1} = r h_k .$$

With $A_k = B_k^{-1}$, the localizing ellipsoid is $\Phi_k = \{x : \|A_k(x - x_k)\| \le (n+1)h_k\}$, and the dimension-dependent ratio is

$$q_n = \sqrt{\frac{n-1}{n+1}}\left(\frac{n}{\sqrt{n^2-1}}\right)^n < 1 .$$

Three problems produce such a field: minimizing a convex $f$ on a ball (the field (3.62), a subgradient inside the ball and the outward radial direction outside); the convex program $\min f_0$ s.t. $f_i \le 0$ (the field (3.65), a subgradient of the objective at feasible points and of a most violated constraint otherwise); and a convex–concave saddle point problem (the field $\{g_f^x, -g_f^y\}$).

Formalization targets

Goal: Theorem 3.14 (p. 86)

For every $k$,

$$\|A_k(x_k - x^*)\| \le h_k (n+1), \tag{3.61}$$

that is, $x^* \in \Phi_k$. The statement holds for every field $g$ satisfying the monotonicity condition at $x^*$; nothing about continuity or convexity is assumed.

Milestones

1. Eq. (3.4): $\|R_\alpha(\xi)x\| = \sqrt{\|x\|^2 + (\alpha^2-1)(x,\xi)^2}$ for unit $\xi$.
2. Volume of $\Phi_k$ (p. 87): $(n+1)h_k = R r^k$ and $v(\Phi_k) = v_0 R^n r^{nk}/\det A_k$, $v_0$ the volume of the unit ball.
3. Volume ratio (p. 87–88): $v(\Phi_{k+1}) = q_n\, v(\Phi_k)$ with $q_n < 1$, the volumes being positive and finite.
4. Eq. (3.62): the ball field satisfies $(g(x), x - x^*) \ge 0$.
5. Eq. (3.65): the convex-programming field satisfies $(g(x), x - x^*) \ge 0$.
6. Saddle point field (p. 90): $(g(z), z - z^*) \ge 0$.

Significance

The result. Theorem 3.14 with the volume identity says that after $k$ steps the solution is confined to an ellipsoid of volume $q_n^k$ times that of the initial ball, for any field of the above kind. Milestones 4–6 turn this into localization guarantees for constrained convex minimization, general convex programming and convex–concave saddle points, with a rate that depends only on the dimension. The same localization underlies the complexity bounds of the ellipsoid method for linear programming and the polynomial equivalence of separation and optimization.

Formalizing it. The results are classical and proved in the book. No machine-checked proof of the space-dilation form of the method is known to exist. The platform already contains a proved version of the Bertsimas–Tsitsiklis form (LinearOptimization.ellipsoid_update_halfspace_subset, LinearOptimization.ellipsoid_update_volume_lt: a half-ellipsoid $E(z, D) \cap \{a^Tx \ge a^Tz\}$ is covered by an updated ellipsoid whose volume is smaller by a factor below $e^{-1/(2(n+1))}$), and Khachiyan's feasibility algorithm (SmaleNinth.khachiyan_ellipsoid_decides). Those statements are about a center/shape-matrix update and give a volume inequality; this mission is about the iterates of Shor's matrix recursion $B_{k+1} = B_k R_\beta(\xi_k)$ and the exact ratio $q_n$. Relating the two parametrizations ($D_k = (n+1)^2 h_k^2 B_k B_k^T$) is a welcome side result.

Difficulty

The obvious approach, tracking the ellipsoid through the center and shape matrix and invoking a minimum-volume covering argument, is not what the algorithm computes: here the iterate is updated through the factor $B_k$ and the stepsize $h_k$ is fixed in advance, independent of the field, so the induction must be carried out in the transformed coordinates $z_k = A_k(x_k - x^*)$ in which the ellipsoid is a ball. The difficulty is that the monotonicity condition gives only the sign of one inner product, $(z_k, \xi_k) \ge 0$, while the norm of $z_{k+1}$ depends on both $(z_k,\xi_k)$ and $\|z_k\|$; the constants $\beta$ and $r$ are exactly those for which the resulting quadratic estimate closes. For the volume identity, the main technical step is computing $\det R_\beta(\xi) = \beta$ and the Lebesgue measure of a linear image of a ball in EuclideanSpace.

Formalization scope

• $E_n$ is EuclideanSpace ℝ (Fin n); matrices act through Matrix.toEuclideanLin; $B_k^*$ is the transpose. $R_\alpha(\xi)$ is the matrix $I + (\alpha-1)\xi\xi^T$ (the book's property 10); for unit $\xi$ this is the operator of the book's definition.
• The algorithm is the definition ellipsoidMethod g R x₀ : ℕ → EllState n, with state $(x_k, B_k, h_k)$, $B_0 = I$, $h_0 = R/(n+1)$. If $g(x_k) = 0$ the state is repeated from then on (the book stops); the normalization in (3.57) is only performed when $g(x_k) \ne 0$. $A_k$ is the matrix inverse of $B_k$, which is nonsingular.
• Theorem 3.14 is stated for $n > 1$, $R > 0$, $x^*$ with $\|x_0 - x^*\| \le R$ and $(g(x), x - x^*) \ge 0$ for all $x$, for every $k$. The book's additional assumption that $g(x) \ne 0$ for $x \ne x^*$ is not used by its proof and is omitted. The iterates are those of the recursion; a statement about an arbitrary ellipsoid containing $x^*$, or about an arbitrary invertible matrix in place of $B_k$, would not be this theorem and is ruled out by the definitions.
• Volumes are Lebesgue measure in ℝ≥0∞. The ellipsoid is given for an arbitrary center (the book writes $x^*$ in one place and $x_k$ in another; the volume is the same). The volume formula requires the first $k$ iterations to have been performed ($g(x_j) \ne 0$, $j < k$); the ratio requires iteration $k+1$ to be performed.
• The printed chain on p. 87 has misprints (exponents $2$ and $1$ on $n/\sqrt{n^2-1}$ where $n$ is meant, and $(n-1)/(n+1)$ for $\sqrt{(n-1)/(n+1)}$); the statement follows the value of $q_n$ given on p. 88. The estimate for $x \notin S(x_0,R)$ before (3.62) has a sign misprint; only the conclusion is stated.
• Needed infrastructure: determinant of a rank-one perturbation of the identity ($\det(I + c\,\xi\xi^T) = 1 + c\|\xi\|^2$, available in Mathlib as the matrix determinant lemma), the measure of a linear image (MeasureTheory.Measure.addHaar_image_linearMap), and elementary real inequalities for $q_n < 1$. The space-dilation lemmas are reusable in the other Shor missions on SDG methods and the $r$-algorithm.

Selected references

• N. Z. Shor, Minimization Methods for Non-Differentiable Functions, Springer Series in Computational Mathematics 3, Springer, 1985, §3.8. https://doi.org/10.1007/978-3-642-82118-9
• D. B. Yudin and A. S. Nemirovskii, Informational complexity and efficient methods for the solution of convex extremal problems, Ekonomika i Matematicheskie Metody 12 (1976), 357–369 (English translation: Matekon 13 (1977), 25–45).
• N. Z. Shor, Cut-off method with space extension in convex programming problems, Cybernetics 13 (1977), 94–96.
• L. G. Khachiyan, A polynomial algorithm in linear programming, Soviet Mathematics Doklady 20 (1979), 191–194.
• R. G. Bland, D. Goldfarb and M. J. Todd, The ellipsoid method: a survey, Operations Research 29 (1981), 1039–1091. https://doi.org/10.1287/opre.29.6.1039
• M. Grötschel, L. Lovász and A. Schrijver, Geometric Algorithms and Combinatorial Optimization, Springer, 1988. https://doi.org/10.1007/978-3-642-97881-4