Convex Optimization

Motivation

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

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

a smooth term $\Phi$ plus a nonsmooth but "simple" term $\Psi$, such as an $\ell^1$ penalty or the indicator function of a convex constraint set. The forward-backward method alternates a gradient step on $\Phi$ with a proximal step on $\Psi$. Beck and Teboulle's FISTA (Beck–Teboulle 2009) combined it with Nesterov's acceleration and improved the worst-case rate for function values from $\mathcal O(k^{-1})$ to $\mathcal O(k^{-2})$. Whether the iterates of the accelerated scheme converge at all, and not just their function values, remained unsettled for a long time; in the words of Attouch and Peypouquet, it "puzzled researchers for over two decades".

Timeline.

• 1967: Opial proves that weak convergence of a sequence in a Hilbert space follows from two facts, the convergence of its distance to every point of a target set and the location of its weak cluster points in that set (Opial 1967).
• 2009: Beck and Teboulle introduce FISTA, with an $\mathcal O(k^{-2})$ rate for function values.
• 2014: Su, Boyd and Candès read the accelerated method as a discretization of the ODE $\ddot x + \frac{\alpha}{t}\dot x + \nabla\Theta(x) = 0$ (Su–Boyd–Candès 2014).
• 2014–2015: for the variant with inertial coefficient $\frac{k-1}{k+\alpha-1}$ and $\alpha > 3$, Chambolle and Dossal (2015) and, independently, Attouch, Chbani, Peypouquet and Redont (arXiv:1507.04782) prove weak convergence of the iterates.
• 2016: Attouch and Peypouquet (arXiv:1510.08740, SIAM J. Optim. 26(3), 2016) prove the $o(k^{-2})$ rate and, as Theorem 3, give a short proof of weak convergence from the same energy estimates.

The case $\alpha = 3$, the original choice of FISTA, is not covered by this result.

Setting

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

• $\Psi : \mathcal H \to \mathbb R \cup \{+\infty\}$ is proper (finite somewhere), lower-semicontinuous and convex. The value $+\infty$ matters: indicator functions of closed convex sets are the main example.
• $\Phi : \mathcal H \to \mathbb R$ is convex and continuously differentiable, and its gradient $\nabla\Phi$ is $L$-Lipschitz continuous.
• $\Theta = \Psi + \Phi$, and $S = \operatorname{argmin}\Theta$ is assumed nonempty.
• For $s > 0$, the proximal map $\operatorname{prox}_{s\Psi}(x)$ is the unique minimizer of $y \mapsto \Psi(y) + \frac{1}{2s}\|y - x\|^2$.

Given $\alpha > 0$ and $s > 0$, algorithm (2) generates $(x_k)$ by

$$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).$$

The auxiliary sequence (6) is $z_k = x_k + \frac{k-1}{\alpha-1}(x_k - x_{k-1})$. For a point $x^*$, the proof of Theorem 3 uses

$$\delta_k = (k-1)\big[\|x_k - x^*\|^2 - \|x_{k-1} - x^*\|^2\big] + (\alpha-1)\|x_k - x^*\|^2 .$$

A sequence $(x_k)$ converges weakly to $\bar x$, written $x_k \rightharpoonup \bar x$, if $\langle x_k, y\rangle \to \langle \bar x, y\rangle$ for every $y \in \mathcal H$.

In Lean these are theta, extrap ($y_k$), IsAccelFBRun and zSeq (shared definitions in the namespace NesterovFB.Rates), and deltaSeq in the namespace NesterovFB.Weak. The predicates IsProperClosedConvex and IsProx and the predicate WeakTendsto are reused from the platform.

Formalization targets

Goal: Theorem 3 (p. 5)

Under the assumptions above, with $\alpha > 3$ and $0 < s < 1/L$,

$$\exists\, \bar x \in S:\quad x_k \rightharpoonup \bar x .$$

The limit is required to be a minimizer of $\Theta$; weak convergence to an arbitrary point would be a weaker statement.

Milestones (proof of Theorem 3, p. 5)

For every $x^* \in S$ and $k \ge 1$:

$$\|x_{k+1} - x^*\|^2 \le \|y_k - x^*\|^2,$$ $$\delta_{k+1} - \delta_k \le 2(k+\alpha-1)\,\|x_k - x_{k-1}\|^2,$$ $$\lim_{k\to\infty}\|z_k - x^*\| \text{ exists},\qquad \lim_{k\to\infty}\|x_k - x^*\| \text{ exists}.$$

Significance

The result. Theorem 3 shows that the accelerated forward-backward method with $\alpha > 3$ behaves like the unaccelerated method in one important respect: its iterates converge to a solution rather than just producing small function values. In infinite-dimensional settings (inverse problems, PDE-constrained optimization, signal recovery in function spaces) weak convergence is the natural notion, and strong convergence can fail. The four milestones are the steps that matter for other analyses too: a "Fejér-type" step inequality from the extrapolated point, and the convergence of the distance to every minimizer.

Formalizing it. The result is proved in the literature; no machine-checked proof of it is known. A complete formalization needs the energy estimates of the paper's §1.1–1.2 (summability of $k\|x_k - x_{k-1}\|^2$, boundedness of $(z_k)$, and the convergence of $k^2\|x_{k+1}-x_k\|^2 + (k+1)^2(\Theta(x_{k+1}) - \min\Theta)$), which a companion mission poses separately, and Opial's lemma, which is not in Mathlib.

Difficulty

The unaccelerated forward-backward method is Fejér monotone: $\|x_{k+1} - x^*\|$ decreases for every minimizer $x^*$, and Opial's lemma applies directly. The accelerated iterates are not Fejér monotone. The first milestone only compares $x_{k+1}$ with the extrapolated point $y_k$, and $\|y_k - x^*\|$ can exceed $\|x_k - x^*\|$ by the inertial term. The quantity $\|x_k - x^*\|^2$ therefore satisfies only a second-order inequality with coefficients that grow in $k$. The obvious attempt, to show that $\|x_k - x^*\|$ is eventually monotone and apply Opial's lemma as for the unaccelerated method, fails.

A second difficulty is the passage from distances to weak convergence: in a Hilbert space this needs the weak sequential compactness of bounded sets and the weak lower-semicontinuity of $\Theta$ (to place weak cluster points in $S$).

Formalization scope

• $\mathcal H$ is a general real Hilbert space (InnerProductSpace ℝ H, CompleteSpace H), not $\mathbb R^n$; in finite dimensions weak and strong convergence coincide and the goal would be a different, weaker theorem.
• $\Psi$ and $\Theta$ take values in EReal; $\Psi$ is never assumed real-valued.
• $L$ is a nonnegative real and "$0 < s < 1/L$" is written $0 < s$, $sL < 1$, which also allows $L = 0$.
• $\operatorname{prox}_{s\Psi}$ is a map $P$ given with its minimization property; under the hypotheses it is unique, so this is the proximal map.
• The run starts at $k = 1$ with $x_0, x_1$ arbitrary. At $k = 1$ the inertial coefficient vanishes, so $x_0$ never matters; every sequence the paper generates satisfies the predicate.
• "The limit exists" means a real limit. Weak convergence is the platform's WeakTendsto, not convergence in norm.
• The milestones keep the standing hypothesis $\alpha > 3$ of Theorem 3, although the first two do not need it.

A trivializing formalization is ruled out by the sanity check: the hypotheses are jointly satisfiable, and the goal asks for weak convergence to a minimizer, so neither a vacuous hypothesis nor an arbitrary limit point is admitted.

Infrastructure that a complete development needs, and that is reusable beyond this mission: the descent inequality (9) of the proximal-gradient operator, Opial's lemma in a real Hilbert space, the weak lower-semicontinuity of proper lower-semicontinuous convex functions, and the lemma that a bounded real sequence whose positive increments are summable converges. Contributions of any of these are welcome.

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. https://arxiv.org/abs/1510.08740 (v4 is the source of this mission)
• H. Attouch, Z. Chbani, J. Peypouquet, P. Redont, Fast convergence of inertial dynamics and algorithms with asymptotic vanishing viscosity, Math. Program. 168, 2018. https://arxiv.org/abs/1507.04782
• A. Chambolle, C. Dossal, On the convergence of the iterates of the "fast iterative shrinkage/thresholding algorithm", J. Optim. Theory Appl. 166, 2015. https://doi.org/10.1007/s10957-015-0746-4
• A. Beck, M. Teboulle, A fast iterative shrinkage-thresholding algorithm for linear inverse problems, SIAM J. Imaging Sci. 2(1):183–202, 2009. https://doi.org/10.1137/080716542
• Z. Opial, Weak convergence of the sequence of successive approximations for nonexpansive mappings, Bull. Amer. Math. Soc. 73:591–597, 1967. https://doi.org/10.1090/S0002-9904-1967-11761-0
• W. Su, S. Boyd, E. J. Candès, A differential equation for modeling Nesterov's accelerated gradient method: theory and insights, NIPS 2014. https://arxiv.org/abs/1503.01243

Motivation

Convex duality in optimization rests on one transformation: to a convex function $f$ one associates the function $\varphi(\xi) = \sup_x(\Sigma x\xi - f(x))$, which records, for each slope $\xi$, the best affine lower bound of $f$ with that slope. Lagrangian duality, the duality theory of linear and conic programming, the analysis of first-order methods through smoothness and strong convexity of conjugates, and the dual representations of risk measures and divergences all read off properties of $f$ from properties of $\varphi$. Each of these uses needs one fact: that the transformation loses no information, so that applying it twice returns $f$.

That fact, for functions on $\mathbb R^n$, is the theorem of W. Fenchel's five-page note On conjugate convex functions (Canad. J. Math. 1 (1949) 73–77). Its timeline:

• 1912. W. H. Young proves the inequality $ab \le F(a) + G(b)$ for a pair of mutually inverse increasing functions $F', G'$ of one variable (Proc. R. Soc. Lond. A 87, 1912).
• 1949. Fenchel defines the conjugate of a convex function on a convex subset of $\mathbb R^n$, without any differentiability, and proves that conjugation is a symmetric correspondence on convex functions that are semi-continuous from below and whose domain is closed relative to the function (Fenchel 1949).
• 1965. J.-J. Moreau develops conjugation for convex functions with values in $(-\infty, +\infty]$ on a real Hilbert space, together with the proximal map (Bull. SMF 93, 1965); the biconjugation theorem in this generality is called the Fenchel–Moreau theorem.
• 1970. R. T. Rockafellar's Convex Analysis makes the conjugate the central object of finite-dimensional convex analysis (Princeton, 1970).

Setting

Points of $\mathbb R^n$ are $x = (x_1,\dots,x_n)$, and $\Sigma x\xi = x_1\xi_1 + \dots + x_n\xi_n$.

A standing pair $(G, f)$ consists of a set $G \subseteq \mathbb R^n$ and a real function $f$ defined in $G$ such that

1. $G$ is nonempty and convex;
2. $f$ is convex on $G$: $f((1-\theta)x' + \theta x'') \le (1-\theta)f(x') + \theta f(x'')$ for $x', x'' \in G$, $0 < \theta < 1$;
3. $f$ is semi-continuous from below on $G$: $\liminf_{x \to x^*,\, x \in G} f(x) \ge f(x^*)$ for $x^* \in G$;
4. $G$ is closed relative to $f$: $f(x) \to +\infty$ as $x \to x^*$ within $G$, for every boundary point $x^*$ of $G$ not in $G$.

$G$ need be neither open, nor closed, nor bounded. The paper writes the lower limit as a "lim" with a bar under it; the milestone texts write it $\lim_{x\to x^*}$, and it always means $\liminf$.

The conjugate pair of $(G, f)$ is

$$\Gamma = \{\xi \in \mathbb R^n : x \mapsto \Sigma x\xi - f(x) \text{ is bounded above on } G\}, \qquad \varphi(\xi) = \sup_{x \in G}\bigl(\Sigma x\xi - f(x)\bigr)\quad (\xi \in \Gamma).$$

The same construction applied to $(\Gamma, \varphi)$ gives the pair $(G^*, f^*)$, with $f^*(x) = \sup_{\xi\in\Gamma}(\Sigma\xi x - \varphi(\xi))$. An interior point of $G$ is a point of the relative interior of $G$, its interior within its affine hull.

In Lean: IsClosedConvexPair G f, conjDomain G f $= \Gamma$, conjFun G f $= \varphi$, and $(G^*, f^*)$ is conjDomain (conjDomain G f) (conjFun G f), conjFun (conjDomain G f) (conjFun G f).

Formalization targets

Goal: Fenchel's theorem (§3, p. 75)

For every standing pair $(G, f)$:

$$(\Gamma,\varphi) \text{ is a standing pair},\qquad \Sigma x\xi \le f(x) + \varphi(\xi)\ \ (x\in G,\ \xi\in\Gamma), \tag{5}$$

with equality for some $\xi \in \Gamma$ at every interior point $x$ of $G$;

$$G^* = G, \qquad f^*(x) = f(x)\ \ (x \in G);$$

and every standing pair $(\Gamma', \varphi')$ whose conjugate pair is $(G, f)$ equals $(\Gamma, \varphi)$.

Milestones, in the order of the proof

1. (5), with no hypothesis on $(G, f)$.
2. $\Gamma \ne \emptyset$, and $\varphi(\xi) = \Sigma x^\circ\xi - f(x^\circ)$ for some $\xi\in\Gamma$ at each interior point $x^\circ$.
3. $\Gamma$ and $\varphi$ are convex.
4. $\varphi$ is semi-continuous from below and $\Gamma$ is closed relative to $\varphi$.
5. (6): $G \subseteq G^*$ and $f^* \le f$ on $G$.
6. Two convex functions, semi-continuous from below on $G$ and equal at the interior points of $G$, are equal on $G$.
7. $f^* = f$ on $G$.
8. (7): $\sup_{\xi\in\Gamma}(\Sigma\xi x^\circ - \varphi(\xi)) = \infty$ for every $x^\circ \notin G$.

Significance

The theorem identifies, among convex functions on convex subsets of $\mathbb R^n$, exactly the class on which conjugation is a bijection and an involution. It is the finite-dimensional base case of the Fenchel–Moreau theorem and underlies Fenchel's duality theorem for $\inf(f - g)$, the conjugate-based optimality conditions of convex programming, and the inversion of gradients of conjugate differentiable convex functions (the paper's §6, the Legendre transformation). The pair form is also how the result is used in practice: the conjugate of a function finite on a set comes with an explicit domain $\Gamma$, and the theorem says that domain determines and is determined by $G$.

The result has been proved for 75 years. What a formalization adds is the theorem in Fenchel's own form: a real-valued function on an explicit convex domain rather than an extended-real function on all of $\mathbb R^n$, the domain identity $G^* = G$ together with the identity of values, and the attainment of equality in (5) at relative-interior points. Machine-checked versions of biconjugation exist in other forms: for extended-real functions on Hilbert spaces, for finite convex functions on all of $\mathbb R^n$, and for functions on a set with a closed restricted epigraph over continuous linear functionals. None of them states the domain identity or the attained equality, and none is in the paper's pair form.

Difficulty

Milestones 1, 3, 4 and 5 follow from the definition of the conjugate alone. The content sits in three places.

• Supporting hyperplanes at relative-interior points. When $G$ is lower-dimensional, the topological interior of $G$ is empty, and a supporting hyperplane must be produced inside the affine hull of $G$ and then extended. Points that are interior to a segment of $G$ but on its relative boundary do not suffice.
• Passing from the interior to the boundary of $G$. Equality $f^* = f$ at interior points does not by itself give equality at boundary points of $G$; it needs semi-continuity from below of both functions and convexity along segments ending at the boundary point.
• $G^* \subseteq G$. Points outside the closure of $G$ and boundary points of $G$ not in $G$ behave differently: a boundary point cannot be separated from $G$ by a hyperplane, and the inclusion there depends on the condition that $G$ be closed relative to $f$. Without that condition the inclusion is false: for $G = (0, 1]$ and $f \equiv 0$, the point $0$ lies in $G^*$.

Formalization scope

• $\mathbb R^n$ is Fin n → ℝ with its product topology, which is the Euclidean one; $\Sigma x\xi$ is Mathlib's x ⬝ᵥ ξ. The paper's $\Sigma\xi x$ in $G^*, f^*$ is ξ ⬝ᵥ x, equal by commutativity.
• $f$ is a total function (Fin n → ℝ) → ℝ, but every hypothesis and conclusion concerns its values on $G$ only: ConvexOn ℝ G f, LowerSemicontinuousOn f G, and Tendsto f (𝓝[G] x) atTop for x ∈ closure G \ G.
• $\varphi$ is the real sSup, which is $0$ on unbounded sets; it is evaluated only on $\Gamma$.
• Three readings are fixed and disclosed in the goal's Formalization Note:
  • (P1) $G$ is nonempty. The paper assumes it tacitly and proves $\Gamma \neq \emptyset$.
  • (P2) Interior points are relative-interior points (intrinsicInterior ℝ G). The paper's segment definition makes the equality clause false, and the topological interior makes it vacuous for lower-dimensional $G$.
  • (P3) Uniqueness is stated as the symmetry gives it. The literal "one and only one $\Gamma$, $\varphi$ with these properties, (5) and equality at interior points" is false: for $G = [0,1]$, $f \equiv 0$, the pair $\Gamma' = \{0\}$, $\varphi'(0) = 0$ also qualifies.
• A statement of the goal that asserts only $f^* = f$ on $G$ and drops $G^* = G$ is a different and much weaker theorem; the goal carries the domain identity.
• Needed infrastructure: supporting hyperplanes to convex sets at relative-interior points (Mathlib has separation theorems for Fin n → ℝ and the intrinsic interior), affine minorants of convex functions on lower-dimensional domains, and the boundary-limit argument of milestone 6. All of these are reusable beyond this mission. Proofs of any milestone, and alternative routes to the goal, are welcome.

Selected references

• W. Fenchel, On conjugate convex functions, Canadian Journal of Mathematics 1 (1949), 73–77. https://doi.org/10.4153/CJM-1949-007-x
• W. H. Young, On classes of summable functions and their Fourier series, Proceedings of the Royal Society of London A 87 (1912), 225–229. https://doi.org/10.1098/rspa.1912.0076
• J.-J. Moreau, Proximité et dualité dans un espace hilbertien, Bulletin de la Société Mathématique de France 93 (1965), 273–299. https://doi.org/10.24033/bsmf.1625
• R. T. Rockafellar, Convex Analysis, Princeton University Press, 1970. https://doi.org/10.1515/9781400873173

Motivation

Many problems in convex optimization, variational inequalities and equilibrium modelling reduce to finding a point $x$ with $0 \in A x + B x$, where $A$ and $B$ are maximal monotone operators on a real Hilbert space $\mathcal H$: for example, minimizing $f + g$ for closed proper convex $f, g$ is the case $A = \partial f$, $B = \partial g$. When the resolvent of $A + B$ is hard to evaluate but the resolvents of $A$ and $B$ separately are easy, one uses a splitting method. Douglas–Rachford splitting, introduced for monotone operators by Lions and Mercier (1979) after an alternating-direction scheme of Douglas and Rachford (1956) for the heat equation, is the most widely used one; through its dual form it underlies the alternating direction method of multipliers (ADMM) used throughout large-scale optimization and statistics.

Eckstein and Bertsekas (MIT report LIDS-P-1919, 1989; Mathematical Programming 55, 1992) showed that Douglas–Rachford splitting is a special case of the proximal point algorithm applied to a single derived operator, the splitting operator $S_{\lambda,A,B}$. This identification lets the convergence theory of the proximal point algorithm transfer to splitting, and yields a generalized method with inexact resolvent evaluations and relaxation.

Timeline.

• Minty (1962): a monotone $T$ is maximal iff $I + T$ is onto.
• Rockafellar (1976): the proximal point algorithm with variable stepsizes and summable errors converges weakly to a zero.
• Lions and Mercier (1979): Douglas–Rachford splitting for maximal monotone $A$, $B$; its map $G_{\lambda,A,B}$ is firmly nonexpansive.
• Gol'shtein and Tret'yakov (1979): relaxed proximal iterations with factors $\rho_k \in (0,2)$, in finite dimension, with a fixed stepsize.
• Eckstein and Bertsekas (1989/1992): the splitting operator; Douglas–Rachford as a proximal point method; the generalized proximal point algorithm and the generalized Douglas–Rachford method, including the case with no solution.

Setting

An operator on $\mathcal H$ is a subset $T \subseteq \mathcal H \times \mathcal H$, with $Tx = \{y \mid (x,y) \in T\}$; it may be multivalued and partially defined. Its domain is $\operatorname{dom} T = \{x \mid Tx \ne \emptyset\}$, its image $\operatorname{im} T$ the projection on the second coordinate, its inverse $T^{-1} = \{(y,x) \mid (x,y) \in T\}$. Scaling and sum are $cT = \{(x, cy)\}$ and $A + B = \{(x, y+z) \mid (x,y) \in A, (x,z) \in B\}$; $I$ is the identity. $T$ is monotone if $\langle x' - x, y' - y\rangle \ge 0$ for all $(x,y),(x',y') \in T$, and maximal monotone if no other monotone operator strictly contains it. The resolvent is $J_{cT} = (I + cT)^{-1}$, and $\operatorname{zer} T = \{x \mid 0 \in Tx\}$. An operator $J$ is firmly nonexpansive if $\|y'-y\|^2 \le \langle x'-x, y'-y\rangle$ for all $(x,y),(x',y') \in J$.

For $\lambda > 0$ the Douglas–Rachford map is $G_{\lambda,A,B} = J_{\lambda A} \circ (2J_{\lambda B} - I) + (I - J_{\lambda B})$, and the splitting operator is

$$S_{\lambda,A,B} = \{(v + \lambda b,\ u - v) \mid (u,b) \in B,\ (v,a) \in A,\ v + \lambda a = u - \lambda b\}.$$

Its zero set is $Z^*_\lambda = \{u + \lambda b \mid b \in Bu,\ -b \in Au\}$.

Formalization targets

Goal: Theorem 7 (generalized Douglas–Rachford splitting)

Let $A$, $B$ be maximal monotone, $\lambda > 0$, and let $\{z^k\}, \{u^k\}, \{v^k\} \subseteq \mathcal H$, $\alpha_k, \beta_k \ge 0$ and $\rho_k$ satisfy

$$\|u^k - J_{\lambda B}(z^k)\| \le \beta_k,\quad \|v^{k+1} - J_{\lambda A}(2u^k - z^k)\| \le \alpha_k,\quad z^{k+1} = z^k + \rho_k (v^{k+1} - u^k),$$

with $\sum \alpha_k < \infty$, $\sum \beta_k < \infty$ and $0 < \inf \rho_k \le \sup \rho_k < 2$. Then

$$\operatorname{zer}(A+B) \ne \emptyset \implies z^k \rightharpoonup z^* \text{ for some } z^* \in Z^*_\lambda,\qquad \operatorname{zer}(A+B) = \emptyset \implies \{z^k\} \text{ unbounded}.$$

Milestones

In the paper's order: Minty's theorem (Theorem 1); properties of firmly nonexpansive operators (Lemma 1); the monotone / firmly nonexpansive correspondence (Theorem 2, Corollaries 2.1–2.3); zeros as fixed points of resolvents (Lemma 2); the generalized proximal point algorithm (Theorem 3): weak convergence to a zero of $T$ under summable errors, relaxation in $(0,2)$ and stepsizes bounded away from $0$, unboundedness when $\operatorname{zer} T = \emptyset$; (maximal) monotonicity of $S_{\lambda,A,B}$ (Theorem 4) and firm nonexpansiveness of its resolvent (Corollary 4.1); $\operatorname{zer} S_{\lambda,A,B} = Z^*_\lambda$ (Theorem 5); and $(I + S_{\lambda,A,B})^{-1} = G_{\lambda,A,B}$ (Theorem 6).

Significance

Theorem 7 gives convergence of Douglas–Rachford splitting with both resolvents evaluated inexactly and with over- or under-relaxation, and it characterizes the case without a solution: the iterates are unbounded exactly when $A + B$ has no zero. The relaxed, inexact form is the one implementations actually run, and through Gabay's identification of ADMM with Douglas–Rachford on the dual it is the basis of the paper's Theorem 8, a convergence theorem for a generalized ADMM. Theorem 3, used to prove Theorem 7, is itself a standard reference form of the inexact relaxed proximal point algorithm.

All results here are proved in the paper (one step in the unbounded case of Theorem 3 rests on results of Rockafellar 1969 and 1970 on sums of maximal monotone operators). As of 2026, neither Douglas–Rachford splitting in this generality nor the generalized proximal point algorithm is formalized in Lean or Mathlib. Mathlib has Hilbert spaces, weak topologies and summability, but no theory of maximal monotone operators, Minty's theorem or resolvents. The mission builds that layer and machine-checks the paper's results on it.

Difficulty

The convergence argument cannot be strong: in infinite dimensions the proximal point algorithm need not converge in norm (Güler 1991), so the conclusion is weak convergence, and identifying the weak limit as a zero requires the weak–strong closedness of the graph of a maximal monotone operator. The maximality halves of Theorems 2 and 4 need Minty's theorem, whose proof requires a nontrivial existence argument (all known proofs use Zorn's lemma or an equivalent). The unbounded case of Theorem 3 is a contradiction argument that truncates $T$ by the subdifferential of the indicator of a ball and invokes two external facts: maximality of the sum of two maximal monotone operators under an interiority condition (Rockafellar 1970), and existence of zeros for maximal monotone operators with bounded domain (Rockafellar 1969). Neither is available in Lean. The natural first idea for Theorem 7, iterating the firm nonexpansiveness of $G_{\lambda,A,B}$, gives neither the error tolerance on both resolvents nor the unbounded case without the full machinery of Theorem 3.

Formalization scope

• $\mathcal H$ is a real inner product space that is complete ([CompleteSpace H]). An operator is a map H → Set H. Monotonicity, maximal monotonicity, $\operatorname{dom}$, $\operatorname{zer}$ and the function-level resolvent predicate IsResolvent are the published definitions ThreeOpSplitting_Convergence_MonotoneOperators; weak convergence is the published WeakTendsto ($\langle z^k, y\rangle \to \langle z^*, y\rangle$ for every $y$).
• §2 notions are graph notions (opResolvent, IsFirmlyNonexpansiveOp, ...), so Theorem 2 and Corollary 2.1 can speak of resolvents that are a priori partial or multivalued. In Theorems 3, 6 and 7 the resolvents are maps $J : \mathcal H \to \mathcal H$ with $\lambda^{-1}(x - J x) \in A(Jx)$ for all $x$, unique by Corollary 2.2.
• $S_{\lambda,A,B}$ is defined by its set formula, not as $G_{\lambda,A,B}^{-1} - I$; with the latter, Theorem 6 and Corollary 4.1 would be unfoldings. Taking free resolvent functions without the IsResolvent hypothesis would make the iteration unrelated to $A$ and $B$; the hypothesis is always present.
• $\inf \rho_k > 0$, $\sup \rho_k < 2$ are encoded as $\exists\, \rho_1, \rho_2$ with $0 < \rho_1 \le \rho_k \le \rho_2 < 2$; $\inf c_k > 0$ as $\exists\, c_0 > 0$, $c_0 \le c_k$. Summability is Summable with nonnegative terms. Sequences start at $k = 0$; $v^0$ is unused. Unboundedness is ¬ Bornology.IsBounded (Set.range z).
• Printed slips corrected and disclosed in the items: Theorem 7 states its sequences in $\mathbb R^n$ (read $\mathcal H$); Theorem 3 prints $(1 - \rho_k) w^k$ (read $\rho_k w^k$, as on p. 9 and in the proof) and $(I + cT)^{-1}$ (read $(I + c_k T)^{-1}$).
• Not included: Corollary 2.4, Corollaries 6.1–6.2 (special cases of Theorem 7), §5 (partial inverses, generalized ADMM). The second sentence of Corollary 6.1 (convergence of $J_{\lambda B}(z^k)$) is deliberately excluded: its argument does not transfer weak convergence (Svaiter 2011).
• Welcome contributions: Minty's theorem in Hilbert space, the resolvent calculus of §2, and weak-limit lemmas (Opial-type arguments) are reusable well beyond this mission.

Selected references

• J. Eckstein and D. P. Bertsekas, On the Douglas–Rachford splitting method and the proximal point algorithm for maximal monotone operators, MIT report LIDS-P-1919, 1989; Mathematical Programming 55 (1992) 293–318. https://doi.org/10.1007/BF01581204
• P.-L. Lions and B. Mercier, Splitting algorithms for the sum of two nonlinear operators, SIAM J. Numer. Anal. 16 (1979) 964–979. https://doi.org/10.1137/0716071
• G. J. Minty, Monotone (nonlinear) operators in Hilbert space, Duke Math. J. 29 (1962) 341–346. https://doi.org/10.1215/S0012-7094-62-02933-2
• R. T. Rockafellar, Monotone operators and the proximal point algorithm, SIAM J. Control Optim. 14 (1976) 877–898. https://doi.org/10.1137/0314056
• O. Güler, On the convergence of the proximal point algorithm for convex minimization, SIAM J. Control Optim. 29 (1991) 403–419. https://doi.org/10.1137/0329022
• B. F. Svaiter, On weak convergence of the Douglas–Rachford method, SIAM J. Control Optim. 49 (2011) 280–287. https://doi.org/10.1137/100788100

Motivation

Many problems in statistics and machine learning minimize an average of $n$ losses, one per data point, plus a regularizer: least squares, logistic regression, and their $\ell_1$- or $\ell_2$-penalized versions. When $n$ is large, a full gradient costs $n$ component gradients, while stochastic gradient descent, which uses one component per step, needs decreasing step sizes and converges slowly. Incremental gradient methods with variance reduction (SAG, SVRG, SDCA, Finito, MISO) use one component gradient per step but converge at the rate of a full-gradient method.

SAGA (Defazio, Bach and Lacoste-Julien, NIPS 2014, arXiv:1407.0202) is a method of this family. It handles a non-smooth regularizer through its proximal operator, and it comes with a guarantee when the losses are convex but not strongly convex. This mission covers that second guarantee, Theorem 2 of the paper. A companion mission covers the linear rate under strong convexity (Theorem 1, Corollary 1).

Timeline.

• 2012: SAG (Le Roux, Schmidt and Bach) gives a linear rate for smooth, strongly convex finite sums. Its analysis does not cover a proximal term.
• 2013: SVRG (Johnson and Zhang) gives a linear rate for the strongly convex case, using periodic full-gradient passes.
• 2013: SDCA (Shalev-Shwartz and Zhang) works on the dual and needs strong convexity.
• 2014: Prox-SVRG (Xiao and Zhang, arXiv:1403.4699) extends SVRG to composite objectives. Its key inequality is reused by SAGA's Theorem 2.
• 2014: SAGA proves both a linear rate under strong convexity and an $O(n/k)$ rate for the averaged iterate under convexity alone, for composite objectives.

Setting

Let $d\ge 0$ and $n\ge 1$. The components $f_1,\dots,f_n:\mathbb R^d\to\mathbb R$ are convex and differentiable, and each gradient $f_i'$ is $L$-Lipschitz ($L>0$). Write

$$f(x)=\frac1n\sum_{i=1}^n f_i(x),\qquad f'(x)=\frac1n\sum_{i=1}^n f_i'(x).$$

The regularizer $h:\mathbb R^d\to\mathbb R$ is convex but possibly non-differentiable. The objective is the composite function $F=f+h$, and $x^*$ is any minimizer of $F$. Minimizers need not be unique, and $f'(x^*)$ need not vanish.

The proximal operator with parameter $\gamma>0$ is

$$\mathrm{prox}_\gamma^h(y)=\arg\min_{x\in\mathbb R^d}\Big\{h(x)+\frac1{2\gamma}\|x-y\|^2\Big\}.$$

SAGA keeps an iterate $x^k$ and a table of points $\phi_1^k,\dots,\phi_n^k$, initialized as $\phi_i^0=x^0$. At step $k+1$ it draws an index $j$ uniformly from $\{1,\dots,n\}$, independently of the past, and sets

$$w^{k+1}=x^k-\gamma\Big[f_j'(x^k)-f_j'(\phi_j^k)+\frac1n\sum_{i=1}^n f_i'(\phi_i^k)\Big],\qquad x^{k+1}=\mathrm{prox}_\gamma^h(w^{k+1}).$$

It then sets $\phi_j^{k+1}=x^k$ and leaves the other table entries unchanged. The averaged iterate is $\bar x^k=\frac1k\sum_{t=1}^k x^t$, which excludes $x^0$.

Formalization targets

Goal: Theorem 2 (p. 11)

With step size $\gamma=1/(3L)$, for every $k\ge1$,

$$\mathbb E\big[F(\bar x^k)\big]-F(x^*)\le\frac{4n}{k}\Big[\frac{2L}{n}\|x^0-x^*\|^2+f(x^0)-\langle f'(x^*),x^0-x^*\rangle-f(x^*)\Big].$$

The expectation is over the indices $j^1,\dots,j^k$. The constants are those printed in the paper.

Milestones (in attack order)

1. Lemma 1 (p. 6) is an inner-product bound for averages of $\mu$-strongly convex functions with $L$-Lipschitz gradients. It is stated for $\mu\ge0$, and Theorem 2 uses the case $\mu=0$.
2. Lemma 2 (p. 7): $\frac1n\sum_i\|f_i'(\phi_i)-f_i'(x^*)\|^2\le 2L\big[\frac1n\sum_i f_i(\phi_i)-f(x^*)-\frac1n\sum_i\langle f_i'(x^*),\phi_i-x^*\rangle\big]$.
3. The bound on $\Delta$ (p. 12). Write $\Delta=-\frac1\gamma(w^{k+1}-x^k)-f'(x^k)$ for the gradient error. For every $\beta>0$, $\mathbb E\|\Delta\|^2\le(1+\beta^{-1})\mathbb E\|f_j'(\phi_j^k)-f_j'(x^*)\|^2+(1+\beta)\mathbb E\|f_j'(x^k)-f_j'(x^*)\|^2$.
4. The prox-SVRG inequality (p. 12): $\alpha\mathbb E\|x^{k+1}-x^*\|^2\le\alpha\|x^k-x^*\|^2-2\alpha\gamma\mathbb E[F(x^{k+1})-F(x^*)]+2\alpha\gamma^2\mathbb E\|\Delta\|^2$.
5. The one-step Lyapunov decrease (p. 12): $\mathbb E[T^{k+1}]-T^k\le-\frac1{4n}\mathbb E[F(x^{k+1})-F(x^*)]$. Here $T(x,\phi)=\frac1n\sum_i f_i(\phi_i)-f(x^*)-\frac1n\sum_i\langle f_i'(x^*),\phi_i-x^*\rangle+(c+\alpha)\|x-x^*\|^2$, with $c=\frac{3L}{2n}$ and $\alpha=\frac{3L}{8n}$.

In milestones 3–5, $\mathbb E$ is the expectation over the single index $j$ of the next step, given the current state.

Significance

The result. Theorem 2 shows that one method, with a step size that depends only on $L$, covers composite problems that are not strongly convex. Examples are $\ell_1$-regularized least squares and logistic regression without a ridge term. On these problems the method converges in expected objective value at rate $O(n/k)$. SAG has no proximal analysis, and SDCA requires strong convexity. With the same step size $1/(3L)$, the paper also states adaptivity to strong convexity, so no strong convexity constant has to be known in advance. The bound is in terms of $T^0$, a quantity computable from the starting point.

Formalizing it. The result is proved on paper, but the proof is not self-contained. Its central inequality (milestone 4) is quoted from the prox-SVRG analysis of Xiao and Zhang, with only the remark that their argument uses $\mathbb E[\Delta]=0$. A machine-checked proof must therefore reconstruct that argument for SAGA's estimator. To our knowledge, no machine-checked proof of SAGA, SVRG or prox-SVRG exists in Lean or Mathlib. The mission also produces reusable statements about convex functions with Lipschitz gradients (Lemmas 1 and 2) and an explicit finite model of a randomized incremental method.

Difficulty

The naive approach applies the non-expansiveness of the proximal operator to $\|x^{k+1}-x^*\|^2$, as in the strongly convex proof. That bounds distances, but it produces no term in $F(x^{k+1})-F(x^*)$. Without strong convexity, the distance terms cannot be traded for function values, so the argument yields no rate.

The function-value term comes from the prox-SVRG inequality (milestone 4), which the paper does not prove. Its difficulty is that $x^{k+1}$ depends on the same random index as $\Delta$, so the cross term between them does not vanish in expectation even though $\mathbb E[\Delta]=0$. A second difficulty is bookkeeping: $w^{k+1}$ uses the old table, the table entry $j$ receives $x^k$ and not $x^{k+1}$, and the constants must make three coefficients vanish exactly. A final step converts the bound on $\frac1k\sum_t\mathbb E[F(x^t)]$ into a bound on $\mathbb E[F(\bar x^k)]$, which requires Jensen's inequality for the convex $F$.

Formalization scope

• Space and indices. Points live in EuclideanSpace ℝ (Fin d). Components are indexed by Fin n (0-based), with $n\ge1$.
• Gradients and smoothness. The gradients are given maps f' with HasGradientAt (f i) (f' i x) x at every point. Smoothness is the Lipschitz bound $\|f_i'(x)-f_i'(y)\|\le L\|x-y\|$.
• Convexity. Convexity is ConvexOn ℝ Set.univ. Lemma 1 uses StrongConvexOn Set.univ μ, whose modulus $\frac\mu2\|x-y\|^2$ is the paper's.
• The regularizer. $h$ is real-valued and convex. Extended-valued regularizers such as indicator functions are outside the statement.
• The proximal map. The proximal operator enters as any map $P$ such that $P(y)$ minimizes $h(z)+\frac1{2\gamma}\|z-y\|^2$ for every $y$. For convex $h$ this determines $P=\mathrm{prox}_\gamma^h$.
• State and expectation. The state is the pair $(x^k,\phi^k)$. The expectation over $k$ steps is the uniform average over the $n^k$ index sequences, which is exactly the law of $k$ independent uniform indices.

Two trivializations are excluded. The averaged-iterate bound carries $k\ge1$, since at $k=0$ the factor $4n/k$ collapses to $0$. The left side is $F$ evaluated at the averaged point, not the average of $F(x^t)$, which is a weaker intermediate step.

A complete development needs the descent lemma and co-coercivity for convex functions with Lipschitz gradients, the characterization and non-expansiveness of the proximal operator, and finite-sum manipulations over index sequences. The lemmas on smooth convex functions and on proximal operators are reusable beyond this mission. Contributions are welcome at every level: proofs of the milestones, a reusable proximal-operator library, and the telescoping argument for the goal.

Selected references

• A. Defazio, F. Bach, S. Lacoste-Julien, SAGA: A Fast Incremental Gradient Method With Support for Non-Strongly Convex Composite Objectives, NIPS 2014. arXiv:1407.0202
• L. Xiao, T. Zhang, A Proximal Stochastic Gradient Method with Progressive Variance Reduction, SIAM J. Optim. 24(4), 2014. arXiv:1403.4699
• N. Le Roux, M. Schmidt, F. Bach, A Stochastic Gradient Method with an Exponential Convergence Rate for Finite Training Sets, NIPS 2012. arXiv:1202.6258
• R. Johnson, T. Zhang, Accelerating Stochastic Gradient Descent using Predictive Variance Reduction, NIPS 2013. NeurIPS proceedings
• S. Shalev-Shwartz, T. Zhang, Stochastic Dual Coordinate Ascent Methods for Regularized Loss Minimization, JMLR 14, 2013. arXiv:1209.1873
• Y. Nesterov, Introductory Lectures on Convex Optimization, Kluwer, 2004. doi:10.1007/978-1-4419-8853-9

Motivation

A Markov decision process (MDP) is solved for one model of its dynamics and rewards, but in practice that model is estimated from data, and a policy that is optimal for the estimate can perform poorly on the true system (Mannor et al., 2007). Robust MDPs address this by evaluating a policy against the worst model in an uncertainty set $\mathcal U$ (Iyengar, 2005; Nilim and El Ghaoui, 2005; Wiesemann, Kuhn and Rustem, 2013). Robust planning, however, solves an inner optimization over $\mathcal U$ at every Bellman update, which is expensive and does not scale to learning settings.

A separate line of work regularizes the policy (entropy, KL, Tsallis penalties) and observes empirically that regularized policies are robust to perturbations (Geist, Scherrer and Pietquin, 2019). Derman, Geist and Mannor (arXiv:2110.06267, NeurIPS 2021) make this precise: for uncertainty sets centred at a nominal model, the robust value function is the solution of a regularized problem posed on the nominal model alone, with a regularizer that is the support function of the uncertainty set. This mission formalizes that equivalence: Proposition 3.1, Theorem 3.1 and Theorem 4.1 of the paper.

Setting

Let $\mathcal S$ and $\mathcal A$ be finite sets of states and actions, $\mathcal A$ nonempty, and $\mathcal X := \mathcal S\times\mathcal A$. Fix a discount factor $\gamma\in(0,1)$ and a strictly positive initial distribution $\mu_0\in\Delta_{\mathcal S}$. A transition kernel $P$ assigns to every pair $(s,a)$ a probability distribution $P(\cdot\mid s,a)$ on $\mathcal S$; a reward is $r\in\mathbb R^{\mathcal X}$. A policy $\pi\in\Delta_{\mathcal A}^{\mathcal S}$ assigns to every state an action distribution $\pi_s$.

For $v\in\mathbb R^{\mathcal S}$ write $r^\pi(s) = \sum_a\pi_s(a)r(s,a)$, $P^\pi(s'\mid s) = \sum_a\pi_s(a)P(s'\mid s,a)$, and define the evaluation Bellman operator

$$T^\pi_{(P,r)}v := r^\pi + \gamma P^\pi v .$$

The inner product on $\mathbb R^{\mathcal S}$ is $\langle v,\mu\rangle = \sum_s v(s)\mu(s)$, and the support function of a set $C\subseteq\mathbb R^{\iota}$ is $\sigma_C(y) = \max_{a\in C}\langle a,y\rangle$.

Given a set $\mathcal U$ of models $(P,r)$, the robust Bellman operator is

$$[T^{\pi,\mathcal U}v](s) := \min_{(P,r)\in\mathcal U}T^\pi_{(P,r)}v(s),$$

and the robust value function $v^{\pi,\mathcal U}$ is its fixed point. Around a nominal model $(P_0,r_0)$, an s-rectangular uncertainty set $\mathcal U = (P_0+\mathcal P)\times(r_0+\mathcal R)$ is given by sets $\mathcal P_s\subseteq\mathbb R^{\mathcal X}$ and $\mathcal R_s\subseteq\mathbb R^{\mathcal A}$, one per state: its models are $P(s'\mid s,a) = P_0(s'\mid s,a)+P_s(s',a)$ and $r(s,a) = r_0(s,a)+r_s(a)$, with $P_s\in\mathcal P_s$ and $r_s\in\mathcal R_s$ chosen independently for each $s$. Finally $[v\cdot\pi_s](s',a) := v(s')\pi_s(a)$.

Formalization targets

Goal: Theorem 4.1 (general robust MDP)

For $\mathcal U = (P_0+\mathcal P)\times(r_0+\mathcal R)$ and every policy $\pi$, $T^{\pi,\mathcal U}$ has a unique fixed point $v^{\pi,\mathcal U}$, and it is the optimal solution of

$$\max_{v\in\mathbb R^{\mathcal S}}\langle v,\mu_0\rangle\quad\text{s.t.}\quad v(s)\le T^\pi_{(P_0,r_0)}v(s)-\sigma_{\mathcal R_s}(-\pi_s)-\sigma_{\mathcal P_s}(-\gamma v\cdot\pi_s)\quad\forall s\in\mathcal S. \tag{2}$$

Milestones

1. Proposition 3.1. For any uncertainty set $\mathcal U = \mathcal P\times\mathcal R$ with $\mathcal P$ a nonempty compact set of kernels and $\mathcal R$ a nonempty compact set of rewards, $v^{\pi,\mathcal U}$ is the optimal solution of the robust program \max_{v}\langle v,\mu_0\rangle\quad\text{s.t.}\quad v\le T^\pi_{(P,r)}v\ \ \forall(P,r)\in\mathcal U. \tag{$P_{\mathcal U}$}
2. Theorem 3.1. For $\mathcal U=\{P_0\}\times(r_0+\mathcal R)$, $v^{\pi,\mathcal U}$ is the optimal solution of $\max_v\langle v,\mu_0\rangle$ s.t. $v(s)\le T^\pi_{(P_0,r_0)}v(s)-\sigma_{\mathcal R_s}(-\pi_s)$ for all $s$.
3. Robust counterpart (proof of Theorem 4.1, App. B.1). For every $v$ and $s$,

$$\max_{(P,r)\in\mathcal U}\{v(s)-r^\pi(s)-\gamma P^\pi v(s)\} = \sigma_{\mathcal P_s}(-\gamma v\cdot\pi_s)+\sigma_{\mathcal R_s}(-\pi_s)+v(s)-T^\pi_{(P_0,r_0)}v(s).$$

Theorem 3.1 is the special case $\mathcal P_s=\{0\}$ of the goal; it is listed separately because it is the paper's statement that policy regularization is equivalent to reward uncertainty.

Significance

The goal says that a robust MDP with s-rectangular uncertainty in both reward and transitions is a regularized MDP on the nominal model, with two regularizers: a policy regularizer $\sigma_{\mathcal R_s}(-\pi_s)$ coming from reward uncertainty, and a regularizer $\sigma_{\mathcal P_s}(-\gamma v\cdot\pi_s)$ coming from transition uncertainty that depends on both the policy and the value. For ball-shaped sets these support functions are explicit ($\alpha^r_s\|\pi_s\|$ and $\alpha^P_s\gamma\|v\|\|\pi_s\|$, Corollary 4.1 of the paper), which leads to the twice regularized (R²) Bellman operators of Section 5 and to robust planning at the cost of non-robust planning. Theorem 3.1 also explains why standard policy regularizers (negative entropy, KL, Tsallis) yield robustness: each is the support function of a reward uncertainty set.

The results are proved in the paper (appendices A.1, A.2, B.1); none has a machine-checked proof. The mission produces formal statements and proofs of the equivalence, the robust Bellman operator's fixed-point theory for stochastic policies and general compact uncertainty sets, and a closed-form robust counterpart that later R² results can import. The paper's printed proof of Proposition 3.1 treats $T^{\pi,\mathcal U}$ as linear in one step; a formal proof settles the statement independently of that step.

Difficulty

The obvious argument reads Proposition 3.1 as linear-programming duality, as for a single MDP. That fails: $T^{\pi,\mathcal U}$ is a minimum of affine maps, hence concave and not affine, and the feasible set of $(P_{\mathcal U})$ is an intersection of infinitely many half-space systems; the argument has to go through monotonicity and contraction of $T^{\pi,\mathcal U}$, which in turn requires every model in $\mathcal U$ to be a genuine transition kernel. For the goal, the paper invokes Fenchel–Rockafellar duality to evaluate the inner maximum; the work in Lean is to separate the maximum over the product set $\mathcal U$ into per-state maxima, which needs the s-rectangular structure and attainment of every maximum (compactness), and to track the index order of the perturbation $P_s(s',a)$ against the kernel $P(s'\mid s,a)$.

Formalization scope

• States and actions are finite types, A nonempty; values are S → ℝ ordered pointwise; a transition array is P : S → A → S → ℝ with P s a s' $=P(s'\mid s,a)$, and the kernel property is the published IsTransitionKernel; $P^\pi(s'\mid s)$ is the published InducedTransition. A policy has π s ∈ stdSimplex ℝ A for every s.
• Perturbations $P_s$ are functions S × A → ℝ indexed $(s',a)$, as in the paper's $\mathbb R^{\mathcal X}$; rewards perturbations are A → ℝ.
• Minima and maxima (in $T^{\pi,\mathcal U}$ and in $\sigma$) are real sInf/sSup. Every theorem assumes the sets nonempty and compact, so these are attained; nothing is quantified over an unbounded set.
• The robust value function is encoded as the fixed point of $T^{\pi,\mathcal U}$, and each theorem asserts its existence and uniqueness. The paper's definition $v^{\pi,\mathcal U}(s)=\min_{(P,r)\in\mathcal U}v^\pi_{(P,r)}(s)$ (p. 4) coincides with it for rectangular sets by a cited result; the proofs use only the fixed-point property. For the non-rectangular sets of Proposition 3.1 the pointwise minimum can be strictly larger than the fixed point and is then not the optimum of $(P_{\mathcal U})$, so the fixed point is the object the proposition is true for.
• "The optimal solution" means: feasible, objective-maximal, and the unique maximizer (uniqueness uses $\mu_0>0$).
• Disclosed hypotheses: $\mathcal U=\mathcal P\times\mathcal R$ with $\mathcal P$, $\mathcal R$ nonempty and compact and every transition in $\mathcal P$ a kernel (Prop. 3.1); $\mathcal P_s$, $\mathcal R_s$ nonempty and compact and every perturbed row $P_0(\cdot\mid s,a)+P_s(\cdot,a)$ in $\Delta_{\mathcal S}$ (Thm 4.1); reward sets rectangular in Thm 3.1, as its proof uses. These are the robust-MDP standing assumptions of p. 4 ($\mathcal P\subseteq\Delta^{\mathcal X}_{\mathcal S}$) and what makes "min" and "max" well defined.
• Not drafted: Corollary 4.1, whose ℓ²-ball $\mathcal P_s$ contains perturbations that leave the simplex, so $P_0+\mathcal P$ is not a set of kernels; Corollary 3.1 and Proposition 3.2 (consequences after the goal; Prop. 3.2 depends on an unspecified policy parametrization).
• A formalization that asserts only that the feasible sets of $(P_{\mathcal U})$ and (2) coincide, or that drops the kernel condition or the existence of the fixed point, does not count: the goal names the robust value function and its optimality.
• "Convex" in the statement of Theorem 4.1 is descriptive and is not part of the formal goal.

Contributions welcome: the monotone-contraction fixed-point lemma for $T^{\pi,\mathcal U}$ and the per-state separation of maxima over rectangular sets are reusable for any robust MDP mission.

Selected references

• E. Derman, M. Geist, S. Mannor, Twice regularized MDPs and the equivalence between robustness and regularization, NeurIPS 2021. arXiv:2110.06267v1
• G. N. Iyengar, Robust dynamic programming, Mathematics of Operations Research 30(2), 2005. doi:10.1287/moor.1040.0129
• A. Nilim, L. El Ghaoui, Robust control of Markov decision processes with uncertain transition matrices, Operations Research 53(5), 2005. doi:10.1287/opre.1050.0216
• W. Wiesemann, D. Kuhn, B. Rustem, Robust Markov decision processes, Mathematics of Operations Research 38(1), 2013. doi:10.1287/moor.1120.0566
• M. Geist, B. Scherrer, O. Pietquin, A theory of regularized Markov decision processes, ICML 2019. PMLR 97
• S. Mannor, D. Simester, P. Sun, J. N. Tsitsiklis, Bias and variance approximation in value function estimates, Management Science 53(2), 2007. doi:10.1287/mnsc.1060.0614

Motivation

Many problems in machine learning and statistics are finite sums: an empirical risk $f(x)=\frac1n\sum_{i=1}^n f_i(x)$ over $n$ data points, often plus a regulariser $h$ such as an $\ell_1$ penalty. When $n$ is large, a full gradient of $f$ costs $n$ component gradients, while stochastic gradient descent uses one component per step but converges only sublinearly because its gradient estimate has non-vanishing variance. Incremental gradient methods with variance reduction keep the per-step cost of one component gradient and still converge linearly on strongly convex problems.

SAGA, introduced by Defazio, Bach and Lacoste-Julien at NIPS 2014 (arXiv:1407.0202), is one of the standard methods of this family, alongside SAG, SVRG, SDCA and Finito/MISO. It keeps a table of past component gradients and handles a non-smooth regulariser through its proximal operator.

Timeline. Le Roux, Schmidt and Bach (2012) gave SAG the first linear rate for strongly convex finite sums at the cost of one gradient per step. Shalev-Shwartz and Zhang (2013) proved linear rates for SDCA, a dual method. Johnson and Zhang (2013) introduced SVRG, with periodic full-gradient passes; Xiao and Zhang (2014) extended it to composite objectives (prox-SVRG). SAGA (2014) combines an unbiased SVRG-style estimator with a SAG-style table, and proves a linear rate in the composite strongly convex case with a simple Lyapunov argument.

Setting

Let $\mathbb R^d$ carry the Euclidean inner product. There are $n\ge1$ differentiable components $f_1,\dots,f_n:\mathbb R^d\to\mathbb R$ with gradients $f_i'$. Each $f_i$ is $\mu$-strongly convex ($\mu>0$): $f_i(ax+by)\le af_i(x)+bf_i(y)-ab\frac\mu2\|x-y\|^2$ for $a,b\ge0$, $a+b=1$. Each gradient is $L$-Lipschitz: $\|f_i'(x)-f_i'(y)\|\le L\|x-y\|$. Write $f=\frac1n\sum_i f_i$ and $f'=\frac1n\sum_i f_i'$. The regulariser $h:\mathbb R^d\to\mathbb R$ is convex, and the goal is to minimise the composite objective $F=f+h$; $x^*$ denotes its minimiser, which is unique.

The proximal operator with step $\gamma>0$ is

$$\operatorname{prox}^h_\gamma(y)=\operatorname*{argmin}_{x}\Big\{h(x)+\tfrac1{2\gamma}\|x-y\|^2\Big\}.$$

SAGA keeps an iterate $x^k$ and points $\phi_1^k,\dots,\phi_n^k$ at which the stored gradients $f_i'(\phi_i^k)$ were taken. It starts from $x^0$ with $\phi_i^0=x^0$. At iteration $k+1$ it draws $j$ uniformly from $\{1,\dots,n\}$, independently of the past, and sets

$$w^{k+1}=x^k-\gamma\Big[f_j'(x^k)-f_j'(\phi_j^k)+\frac1n\sum_{i=1}^n f_i'(\phi_i^k)\Big],\qquad x^{k+1}=\operatorname{prox}^h_\gamma(w^{k+1}),$$

then $\phi_j^{k+1}=x^k$, with every other entry unchanged.

The analysis uses the Lyapunov function

$$T(x,\{\phi_i\})=\frac1n\sum_i f_i(\phi_i)-f(x^*)-\frac1n\sum_i\langle f_i'(x^*),\phi_i-x^*\rangle+c\|x-x^*\|^2 .$$

Formalization targets

Goal: Corollary 1 (p. 8)

With $\gamma=\frac1{2(\mu n+L)}$, for every $k\ge0$,

$$\mathbb E\|x^k-x^*\|^2\le\Big(1-\frac{\mu}{2(\mu n+L)}\Big)^k\Big[\|x^0-x^*\|^2+\frac{n}{\mu n+L}\big(f(x^0)-\langle f'(x^*),x^0-x^*\rangle-f(x^*)\big)\Big],$$

where the expectation is over the indices drawn in the first $k$ iterations. The constants are the paper's.

Theorem 1 (p. 7)

With $\gamma$ as above, $c=\frac1{2\gamma(1-\gamma\mu)n}$ and $\kappa=\frac1{\gamma\mu}$, for every state $(x^k,\{\phi^k_i\})$,

$$\mathbb E\big[T^{k+1}\big]\le\Big(1-\frac1\kappa\Big)T^k ,$$

with the expectation over the next index only.

Supporting lemmas

Lemma 4 (p. 10), a lower bound combining strong convexity and smoothness; Lemma 1 (pp. 6–7), its average over the components; Lemma 2 (p. 7), which bounds the stale-gradient variance by the table part of $T$; and Lemma 3 (p. 7), a second-moment bound for the SAGA step.

Significance

The result. Corollary 1 gives an $\varepsilon$-accurate iterate in expectation after $O\big((n+L/\mu)\log(1/\varepsilon)\big)$ component-gradient evaluations. This is the complexity of full-gradient descent with the condition number decoupled from $n$, and it holds in the composite setting, so it covers the lasso and elastic-net problems that SAG's analysis does not reach. The paper notes that the rate improves on the published rates of SAG and SVRG and is within a factor 2 of SDCA's. Theorem 1 is the template of later Lyapunov analyses of variance-reduced methods.

Formalizing it. The result has been proved since 2014, and no machine-checked proof is known to this mission. The work left is to formalize the known proof: the convexity inequalities (Lemmas 4, 1, 2), the variance computation (Lemma 3), the one-step contraction (Theorem 1), and the passage from conditional to total expectation along the random index sequence (Corollary 1). The paper's Lemma 3 has a sign misprint, which the formalization corrects; see the scope section.

Difficulty

The obvious argument for SGD-type methods bounds $\mathbb E\|x^{k+1}-x^*\|^2$ in terms of $\|x^k-x^*\|^2$ alone. That fails here: the variance of the SAGA estimator depends on the stale table points $\phi_i^k$, which can be far from $x^*$ even when $x^k$ is close. One needs a potential that also measures the table. Balancing the terms of $T$ then requires the four round-bracket coefficients in the paper's display (10) to be non-positive for the specific $\gamma$, $c$ and an auxiliary $\beta=(2\mu n+L)/L$. Checking these coefficients is routine but long algebra in $\mu$, $L$, $n$. The composite case adds one step: since $f'(x^*)\neq0$ in general, the argument goes through the fixed-point identity $x^*=\operatorname{prox}^h_\gamma(x^*-\gamma f'(x^*))$ and the non-expansiveness of the proximal operator, neither of which is a numbered result of the paper.

Formalization scope

• Space and data. The space is EuclideanSpace ℝ (Fin d). Components are indexed by Fin n (0-based), with 0 < n. The gradients $f_i'$ are given maps with HasGradientAt (f i) (f' i x) x. Strong convexity is Mathlib's StrongConvexOn Set.univ μ, whose modulus $\frac\mu2\|x-y\|^2$ is the paper's.

• Regulariser and minimiser. $h$ is real-valued and convex; extended-valued regularisers are out of scope, as on the page. A minimiser $x^*$ of $f+h$ is a hypothesis.

• Proximal operator. It is any map $P$ such that $P(y)$ minimises $h(x)+\frac1{2\gamma}\|x-y\|^2$ for every $y$ (IsProxPoint). The minimiser is unique, so $P$ is $\operatorname{prox}^h_\gamma$.

• State and expectation. The state is the pair $(x,\phi)$. The run after $k$ steps is a deterministic function of the index sequence in Fin k → Fin n. The expectation in Corollary 1 is the average over all $n^k$ sequences, which is exactly the law of $k$ independent uniform indices; no measure theory is involved. Theorem 1's conditional expectation is the average over the next index.

• Constants and corrections. Constants are as printed and fixed, not "for some constant" and not "for all small enough steps". Lemma 4 carries the hypothesis $\mu<L$, which its fractions $1/(L-\mu)$ require. Lemma 3 is stated with $+\gamma f'(x^*)$, as in its proof and its use in Theorem 1; the printed $-\gamma f'(x^*)$ is false whenever $f'(x^*)\neq0$.

• Trivializing formalizations, ruled out. Taking the proximal step as merely non-expansive, fixing an index sequence instead of averaging over all of them, measuring $x^*$ against $f$ instead of $f+h$, or restricting Theorem 1 to reachable states changes the theorem and is excluded.

• Infrastructure. A complete development needs:

  • the co-coercivity inequality for convex functions with Lipschitz gradient;
  • existence, uniqueness and non-expansiveness of the proximal map of a finite convex function;
  • the optimality condition $x^*=\operatorname{prox}^h_\gamma(x^*-\gamma f'(x^*))$;
  • finite-sum variance identities.

  These pieces are reusable well beyond SAGA, by SVRG, SAG and proximal-gradient analyses. Contributions of any of them, or of proofs of the individual milestones, are welcome.

Selected references

• A. Defazio, F. Bach, S. Lacoste-Julien, SAGA: A Fast Incremental Gradient Method With Support for Non-Strongly Convex Composite Objectives, NIPS 2014. arXiv:1407.0202
• N. Le Roux, M. Schmidt, F. Bach, A Stochastic Gradient Method with an Exponential Convergence Rate for Finite Training Sets, NIPS 2012. arXiv:1202.6258
• R. Johnson, T. Zhang, Accelerating Stochastic Gradient Descent using Predictive Variance Reduction, NIPS 2013. NeurIPS proceedings
• L. Xiao, T. Zhang, A Proximal Stochastic Gradient Method with Progressive Variance Reduction, SIAM J. Optim. 24(4), 2014. arXiv:1403.4699
• S. Shalev-Shwartz, T. Zhang, Stochastic Dual Coordinate Ascent Methods for Regularized Loss Minimization, JMLR 14, 2013. arXiv:1209.1873
• Y. Nesterov, Introductory Lectures on Convex Optimization, Kluwer, 2004. doi:10.1007/978-1-4419-8853-9