Convex Optimization

Motivation

Singular value thresholding (SVT) is a first-order method introduced by Cai, Candès and Shen (SIAM J. Optim. 20 (2010)) for recovering a low-rank matrix from incomplete or indirect information. Its basic form, for matrix completion, alternates a soft-thresholding of singular values with a gradient step on a dual variable, and needs only one sparse singular value decomposition per iteration. That is what made nuclear-norm heuristics usable on matrices with tens of thousands of rows and columns, where interior-point methods for the equivalent semidefinite program do not fit in memory.

Matrix completion is only one constraint set. In applications the data are noisy linear measurements $b = \mathcal A(M) + z$, and the constraint takes the form of componentwise error bounds or norm balls around the data (§3.3 of the paper). Section 3.2 of the paper extends the method to a general finite family of convex constraints, and §4.2 proves that the extended iteration converges. This mission formalizes that extension and its convergence theorem, Theorem 4.4.

Setting

Let $n_1, n_2, m$ be natural numbers and $\mathbb R^{n_1\times n_2}$ the real $n_1\times n_2$ matrices, with the Frobenius inner product $\langle X, Y\rangle = \sum_{i,j} X_{ij}Y_{ij}$ and norm $\|X\|_F = \sqrt{\langle X, X\rangle}$. The nuclear norm $\|X\|_*$ is the sum of the singular values of $X$. For a fixed $\tau > 0$ the objective is

$$f_\tau(X) = \tau\|X\|_* + \tfrac12\|X\|_F^2 .$$

A matrix $Z$ is a subgradient of a function $g$ at $X_0$, written $Z\in\partial g(X_0)$, if $g(X)\ge g(X_0) + \langle Z, X - X_0\rangle$ for all $X$.

Let $f_1,\dots,f_m:\mathbb R^{n_1\times n_2}\to\mathbb R$ be convex and put $\mathcal F(X) = (f_1(X),\dots,f_m(X))\in\mathbb R^m$. On $\mathbb R^m$, $\langle u, v\rangle = \sum_i u_iv_i$ and $\|v\|$ is the Euclidean norm. The constrained problem is

$$\text{(3.4)}\qquad \text{minimize } f_\tau(X)\quad\text{subject to } f_i(X)\le 0,\ i=1,\dots,m,$$

with Lagrangian $\mathcal L(X, y) = f_\tau(X) + \langle y, \mathcal F(X)\rangle$ for $y\ge 0$. A pair $(X^\star, y^\star)$ with $y^\star\ge0$ is primal-dual optimal if it is a saddle point:

$$\mathcal L(X^\star, y)\le \mathcal L(X^\star, y^\star)\le \mathcal L(X, y^\star)\qquad\text{for all } y\ge 0,\ X .$$

The paper's standing assumption "strong duality holds" is the existence of such a pair.

The iteration (3.5) starts from $y^0 = 0$ and, for step sizes $\delta_k$, sets for $k = 1, 2, \dots$

$$X^k = \arg\min_X\{f_\tau(X) + \langle y^{k-1}, \mathcal F(X)\rangle\},\qquad y^k = [\,y^{k-1} + \delta_k\mathcal F(X^k)\,]_+ ,$$

where $x_+$ has entries $\max(x_i, 0)$. It is Uzawa's method for (3.4): an exact minimization in the primal variable followed by a projected ascent step on the dual. When $\mathcal F(X) = b - \mathcal A(X)$ is affine, the minimization is a singular value thresholding step, which gives the algorithm its name.

The analysis of §4.2 assumes $\mathcal F$ is Lipschitz in the sense

$$\text{(4.2)}\qquad \|\mathcal F(X) - \mathcal F(Y)\|\le L\,\|X - Y\|_F\quad\text{for all } X, Y,$$

for a constant $L\ge 0$.

Formalization targets

Goal: Theorem 4.4 (p. 1969)

If $0 < \inf_k\delta_k\le\sup_k\delta_k < 2/L^2$ and strong duality holds, then the sequence $X^k$ of (3.5) converges to the unique solution of (3.4):

$$\exists!\,X^\star\ \text{solving (3.4)},\qquad \lim_{k\to\infty} X^k = X^\star .$$

Milestones, in the order the proof uses them

• Lemma 4.1 (p. 1968): $\langle Z - Z', X - X'\rangle\ge\|X - X'\|_F^2$ for $Z\in\partial f_\tau(X)$, $Z'\in\partial f_\tau(X')$.
• Lemma 4.3 (p. 1969): for a primal-dual optimal pair and each $\delta > 0$, $y^\star = [y^\star + \delta\mathcal F(X^\star)]_+$.
• Eq. (4.4) (p. 1969): there are $Z^k\in\partial f_\tau(X^k)$ and $Z^\star\in\partial f_\tau(X^\star)$ with $\langle Z^k, X - X^k\rangle + \langle y^{k-1}, \mathcal F(X) - \mathcal F(X^k)\rangle\ge 0$ and $\langle Z^\star, X - X^\star\rangle + \langle y^\star, \mathcal F(X) - \mathcal F(X^\star)\rangle\ge 0$ for all $X$.
• Eq. (4.5) (p. 1969): $\langle y^{k-1} - y^\star, \mathcal F(X^k) - \mathcal F(X^\star)\rangle\le -\|X^k - X^\star\|_F^2$.
• Contraction step (p. 1969): $\|y^k - y^\star\|\le\|y^{k-1} - y^\star + \delta_k(\mathcal F(X^k) - \mathcal F(X^\star))\|$.
• Eq. (4.6) (p. 1970): if $2\delta_k - \delta_k^2L^2\ge\beta > 0$ for $k\ge1$, then $\|y^k - y^\star\|^2\le\|y^{k-1} - y^\star\|^2 - \beta\|X^k - X^\star\|_F^2$.

Significance

The result. Theorem 4.4 is the convergence guarantee for SVT beyond matrix completion. The componentwise error bounds of (3.8), whose SVT iteration is (3.9), are finitely many affine constraints and fall under it directly, as does any finite family of Lipschitz convex constraints, for instance a Frobenius-norm ball around the data. The conic variants of §3.3 ((3.11)–(3.13)) project the dual variable onto a cone rather than onto the nonnegative orthant and are not covered by the theorem as stated. Together with Theorem 3.1 of the same paper, which says that the solution of (3.4) tends to the minimum-nuclear-norm solution as $\tau\to\infty$, it justifies using SVT as a solver for nuclear-norm minimization under general convex constraints.

Formalizing it. The theorem is proved in the paper, with two steps delegated to the literature: Lemma 4.3 cites [31], and the concluding step reads "the conclusion is as before". Its proof is short but relies on convex-analytic facts that are standard on paper and missing, in this form, from Mathlib: subgradients of the nuclear norm, the subdifferential sum rule for finite convex functions, and nonexpansiveness of the projection onto the nonnegative orthant. No machine-checked proof of this theorem or of Uzawa-type convergence for nuclear-norm objectives is known to exist. The mission produces a complete, checked version of the argument, including the omitted closing step.

Difficulty

The obvious approach is to view (3.5) as projected gradient ascent on the dual function $g(y) = \min_X\mathcal L(X, y)$ and quote the standard convergence theorem for gradient methods with Lipschitz gradients. That does not apply directly: for general convex $f_i$ the dual function need not be differentiable, $\mathcal F(X^k)$ is only a supergradient, and the Lipschitz hypothesis (4.2) is on $\mathcal F$, not on a dual gradient. The proof instead works with the primal-dual pair: it needs first-order optimality conditions (4.4), which require a subdifferential sum rule for $f_\tau + \sum_i y_i f_i$ with nonsmooth $f_i$, and it needs the strong monotonicity of $\partial f_\tau$ (Lemma 4.1), which depends on the description of subgradients of the nuclear norm. A second subtlety is that the theorem asserts convergence of the whole primal sequence to the unique solution, not to some solution along a subsequence, while nothing is claimed about convergence of the dual sequence.

Formalization scope

Matrices are Matrix (Fin n₁) (Fin n₂) ℝ, vectors in $\mathbb R^m$ are Fin m → ℝ, and convergence of matrices is in Mathlib's product topology, which coincides with the Frobenius topology. The nuclear norm is the sum of Mathlib's LinearMap.singularValues of the matrix viewed as a map between Euclidean spaces. Each $f_i$ is a real-valued function with ConvexOn ℝ Set.univ. The iteration is a predicate on sequences indexed by ℕ: the paper's step $k$ produces X (k+1) and y (k+1) from y k with step size δ (k+1), and y 0 = 0. $X^k$ is required to minimize $\mathcal L(\cdot, y^{k-1})$; for $\tau>0$ and convex $f_i$ this minimizer exists and is unique, so the predicate is satisfiable and determines the sequence. The step-size condition is stated as $a\le\delta_k\le C$ for $k\ge1$ with $a>0$ and $C L^2 < 2$, which avoids the division $2/L^2$ (evaluated as $0$ in Lean when $L=0$); for $L=0$ it requires only bounded steps, matching the convention $2/0 = \infty$. Strong duality is the hypothesis that a saddle point exists; Slater's condition is not assumed. The paper's standing assumptions ($\tau>0$, convex $f_i$, and (4.2) where $L$ enters) appear as explicit hypotheses in every statement.

A formalization that assumes convergence or boundedness of the dual iterates, replaces the primal minimization by a closed-form thresholding step (valid only for affine $\mathcal F$), or states only subsequential convergence would not be this theorem; each of these is excluded by the statements above.

A complete development needs: subgradients of the nuclear norm and strong monotonicity of $\partial f_\tau$; existence and characterization of minimizers of strongly convex continuous functions on a finite-dimensional space; the subdifferential sum rule for finite convex functions; complementary slackness from the saddle-point inequalities; and nonexpansiveness of the entrywise positive part. These are reusable beyond this mission, especially for other Uzawa and augmented Lagrangian analyses. Contributions of any of these pieces as separate lemmas are welcome.

Selected references

• J.-F. Cai, E. J. Candès, Z. Shen, A Singular Value Thresholding Algorithm for Matrix Completion, SIAM J. Optim. 20(4):1956–1982, 2010. https://doi.org/10.1137/080738970
• E. J. Candès, B. Recht, Exact Matrix Completion via Convex Optimization, Found. Comput. Math. 9:717–772, 2009. https://doi.org/10.1007/s10208-009-9045-5
• K. J. Arrow, L. Hurwicz, H. Uzawa, Studies in Linear and Nonlinear Programming, Stanford University Press, 1958.
• S. Boyd, L. Vandenberghe, Convex Optimization, Cambridge University Press, 2004. https://doi.org/10.1017/CBO9780511804441

Motivation

Matrix completion asks to recover an $n_1\times n_2$ matrix $M$ from a subset $\Omega$ of its entries. When $M$ has low rank, a standard convex surrogate is to minimize the nuclear norm $\|X\|_*$ (the sum of the singular values) subject to agreeing with $M$ on $\Omega$; Candès and Recht showed that this recovers $M$ exactly under incoherence and sampling conditions (Candès–Recht 2009). Generic interior-point solvers for this semidefinite program do not scale beyond matrices of a few hundred rows.

Cai, Candès and Shen (SIAM J. Optim. 2010) proposed the singular value thresholding (SVT) algorithm: a first-order iteration whose only nonlinear step is a soft-thresholding of singular values, and whose other iterate is a sparse matrix supported on $\Omega$. The algorithm has become a standard baseline in low-rank matrix recovery and a model example of dual (Uzawa-type) methods for nuclear-norm problems. This mission formalizes its convergence theorem.

Setting

All matrices are real. For $X,Y\in\mathbb R^{n_1\times n_2}$ write $\langle X,Y\rangle=\operatorname{trace}(X^*Y)=\sum_{i,j}X_{ij}Y_{ij}$ and $\|X\|_F^2=\langle X,X\rangle$. The nuclear norm $\|X\|_*$ is the sum of the singular values of $X$.

For an index set $\Omega$, the sampling projector $P_\Omega$ keeps the entries with indices in $\Omega$ and sets the others to zero.

A reduced singular value decomposition of a matrix $Y$ of rank $r$ is $Y=U\Sigma V^*$ with $U$ ($n_1\times r$) and $V$ ($n_2\times r$) having orthonormal columns and $\Sigma=\operatorname{diag}(\sigma_1,\dots,\sigma_r)$ with $\sigma_i>0$. For $\tau\ge0$ the singular value shrinkage operator is

$$\mathcal D_\tau(Y)=U\operatorname{diag}\big((\sigma_i-\tau)_+\big)V^*,\qquad t_+=\max(0,t).$$

Fix $\tau>0$, a sequence of step sizes $\{\delta_k\}_{k\ge1}$ and data $M$. The SVT iteration (2.7) starts from $Y^0=0$ and sets, for $k=1,2,\dots$,

$$X^k=\mathcal D_\tau(Y^{k-1}),\qquad Y^k=Y^{k-1}+\delta_k P_\Omega(M-X^k).$$

The proximal problem (2.8) is

$$\text{minimize}\ \ f_\tau(X)=\tau\|X\|_*+\tfrac12\|X\|_F^2\quad\text{subject to}\ \ P_\Omega(X)=P_\Omega(M).$$

More generally, for a linear map $\mathcal A:\mathbb R^{n_1\times n_2}\to\mathbb R^m$ with adjoint $\mathcal A^*$ and spectral norm $\|\mathcal A\|=\sup\{\|\mathcal A(X)\|_{\ell_2}:\|X\|_F=1\}$, and $b\in\mathbb R^m$, problem (3.1) is to minimize $f_\tau(X)$ subject to $\mathcal A(X)=b$, and Uzawa's iteration (3.3) starts from $y^0=0$ and sets $X^k=\mathcal D_\tau(\mathcal A^*(y^{k-1}))$, $y^k=y^{k-1}+\delta_k(b-\mathcal A(X^k))$.

Formalization targets

Goal: Theorem 4.2, second sentence (p. 1968)

If $0<\inf_k\delta_k\le\sup_k\delta_k<2$, then (2.8) has a unique solution $X^\star$ and the SVT iterates satisfy

$$\lim_{k\to\infty}X^k=X^\star .$$

Theorem 4.2, first sentence (p. 1968)

If (3.1) is feasible and $0<\inf_k\delta_k\le\sup_k\delta_k<2/\|\mathcal A\|^2$, then (3.1) has a unique solution and the iterates $X^k$ of (3.3) converge to it.

Supporting results (milestones, in attack order)

1. Well-definedness of $\mathcal D_\tau$ (§2.1, p. 1960): the output does not depend on the chosen SVD.
2. Theorem 2.1 (p. 1960): $\mathcal D_\tau(Y)=\arg\min_X \tfrac12\|X-Y\|_F^2+\tau\|X\|_*$.
3. Sparsity of the iterates (§2.2, p. 1961): since $Y^0=0$, every $Y^k$ vanishes outside $\Omega$.
4. Eq. (2.14) (p. 1964): the minimizers of the Lagrangian $f_\tau(X)+\langle Y,P_\Omega(M-X)\rangle$ are those of $\tau\|X\|_*+\tfrac12\|X-P_\Omega Y\|_F^2$.
5. Lemma 4.1 (p. 1968): for $Z\in\partial f_\tau(X)$, $Z'\in\partial f_\tau(X')$, $\langle Z-Z',X-X'\rangle\ge\|X-X'\|_F^2$.
6. The §3.1 reduction (p. 1964): for a sampling operator, $\mathcal A^*\mathcal A=P_\Omega$ and (3.3) becomes (2.7) under $Y^k=\mathcal A^*(y^k)$.
7. Theorem 4.2, first sentence, as above.

Significance

The theorem certifies that SVT, run with any step sizes in a fixed interval $(0,2)$, computes the unique minimizer of the strongly convex surrogate (2.8). Together with the separate fact that the solution of (2.8) tends to the minimum nuclear norm completion as $\tau\to\infty$ (the paper's Theorem 3.1, a companion mission), this is what justifies using SVT as a solver for nuclear-norm matrix completion. Theorem 2.1, the proximal characterization of singular value soft-thresholding, is used throughout the literature on proximal methods for low-rank problems.

The paper's proof of Theorem 4.2 consists of the reduction to Uzawa's method and a citation of a general convergence theorem for projected gradient methods on the dual. The formalization produces a self-contained, machine-checked chain: the proximal characterization of $\mathcal D_\tau$, the Lagrangian identity, strong monotonicity of $\partial f_\tau$, and the convergence argument itself. To our knowledge none of these results has a machine-checked proof; Mathlib at the pinned revision has singular values of linear maps but no SVD structure, no nuclear norm and no subgradient calculus.

Difficulty

Nothing in the iteration is a gradient step of a smooth function in $X$: the $X$-update is a nonsmooth proximal map, and the convergence of $X^k$ is not visible from the recursion itself. The paper's argument cites a general theorem on projected gradient methods ([25, Theorem 2.1]) and takes for granted that "strong duality holds" for (2.8) (p. 1963), so the existence of a Lagrange multiplier is part of what must be formalized. Theorem 2.1 depends on the subdifferential of the nuclear norm, which Mathlib does not provide, and therefore on the singular value decomposition and the duality between the nuclear and spectral norms. Convergence of objective values or of a subsequence would not suffice: the target is convergence of the whole sequence $X^k$ to the unique solution.

Formalization scope

Matrices are Matrix (Fin n₁) (Fin n₂) ℝ; convergence is Mathlib's topology on matrices, which coincides with the Frobenius-norm topology. $\langle\cdot,\cdot\rangle$ and $\|\cdot\|_F$ are defined entrywise; $\|X\|_*$ is the sum of Mathlib's LinearMap.singularValues of $X$ viewed as a map $\mathbb R^{n_2}\to\mathbb R^{n_1}$. The shrinkage operator is a relation IsShrink τ Y X defined, as in (2.1)–(2.2), through some reduced SVD of $Y$; well-definedness is a milestone. It is not defined as the minimizer of (2.3), which would make Theorem 2.1 definitional. Linear maps $\mathcal A$ are given by matrices $A_1,\dots,A_m$ with $\mathcal A(X)_i=\langle A_i,X\rangle$, and a sampling operator by an injective enumeration of $\Omega$. Subgradients are those of (2.4).

Sequences are indexed by $\mathbb N$: Lean's step $k+1$ is the paper's step $k$, so $X^0$ and $\delta_0$ are unused. Committed conventions:

• $Y^0=0$ and $y^0=0$ are hypotheses; with a start that is nonzero outside $\Omega$ the iterates converge to a different matrix.
• The standing $\tau>0$ is kept, except in Theorem 2.1 and the well-definedness statement, which are printed for $\tau\ge0$.
• The step-size conditions are explicit bounds $a>0$, $C$ with $a\le\delta_k\le C$ for $k\ge1$, together with $C<2$, respectively $C\|\mathcal A\|^2<2$. The multiplicative form avoids Lean's $x/0=0$: for $\mathcal A=0$ the condition does not become unsatisfiable.
• Feasibility of (3.1) is an added hypothesis of Theorem 4.2's first sentence, since "the unique solution" presupposes it.
• "Converges to the unique solution" is stated as existence and uniqueness of the solution together with convergence of the whole sequence to it.

A small unprinted helper, $\|\mathcal A\|\le1$ for sampling operators, is included to pass from the first sentence of Theorem 4.2 to the second; it is not a milestone. The subdifferential formula (2.6) of the nuclear norm and the Fejér-type condition of §5.1.2 are not stated. Reusable infrastructure welcome from solvers: existence and uniqueness properties of the reduced SVD, the nuclear/spectral norm duality, the subdifferential of the nuclear norm, and a general convergence theorem for Uzawa's method with a strongly convex objective.

Selected references

• J.-F. Cai, E. J. Candès, Z. Shen, A Singular Value Thresholding Algorithm for Matrix Completion, SIAM J. Optim. 20(4):1956–1982, 2010. https://doi.org/10.1137/080738970
• E. J. Candès, B. Recht, Exact Matrix Completion via Convex Optimization, Found. Comput. Math. 9:717–772, 2009. https://doi.org/10.1007/s10208-009-9045-5
• K. J. Arrow, L. Hurwicz, H. Uzawa, Studies in Linear and Non-Linear Programming, Stanford University Press, 1958.

Motivation

Many problems in signal processing, statistics and machine learning minimise a sum of three convex terms: a smooth data-fit term and two nonsmooth regularisers or constraints, each of which is easy to handle on its own (through its proximal map) but not in combination. Examples are constrained sparse regression, matrix completion with a nuclear-norm penalty and box constraints, and support-vector machines with a norm penalty. Davis and Yin (Set-Valued Var. Anal. 25 (2017)) introduced a three-operator splitting scheme that evaluates each proximal map and the gradient of the smooth term once per iteration and reduces to Douglas–Rachford splitting (Lions and Mercier 1979) and forward–backward splitting as special cases. Section 3 of that paper gives the objective-error rates of the scheme on convex problems. This mission formalizes those rates for general convex problems.

Setting

Let $H$ be a real Hilbert space. The problem is

$$\min_{x \in H}\; f(x) + g(x) + h(x), \tag{3.1}$$

where $f, g : H \to (-\infty, +\infty]$ are closed, proper, convex functions (lower semicontinuous, never $-\infty$, finite somewhere, with convex epigraph) and $h : H \to \mathbb R$ is convex and differentiable with $\beta^{-1}$-Lipschitz gradient $\nabla h$, $\beta > 0$.

For $\gamma > 0$ the proximal map $\operatorname{prox}_{\gamma f}(x)$ is the unique minimiser of $y \mapsto f(y) + \frac{1}{2\gamma}\|y - x\|^2$. Algorithm 2 of the paper picks $z^0 \in H$ and $\gamma \in (0, 2\beta)$ and iterates, with relaxation $\lambda_k \equiv 1$,

$$x^k_g = \operatorname{prox}_{\gamma g}(z^k),\qquad x^k_f = \operatorname{prox}_{\gamma f}\big(2x^k_g - z^k - \gamma\nabla h(x^k_g)\big),\qquad z^{k+1} = z^k + x^k_f - x^k_g .$$

Equivalently $z^{k+1} = T z^k$ for the three-operator map

$$T z = \operatorname{prox}_{\gamma f}\big(2\operatorname{prox}_{\gamma g}(z) - z - \gamma\nabla h(\operatorname{prox}_{\gamma g}(z))\big) + z - \operatorname{prox}_{\gamma g}(z).$$

If $z^*$ is a fixed point of $T$, then $x^* = \operatorname{prox}_{\gamma g}(z^*)$ minimises (3.1). The weighted ergodic iterate is

$$\bar x^k_g = \frac{2}{(k+1)(k+2)}\sum_{i=0}^{k} (i+1)\,x^i_g ,$$

and $\bar x^k_f$ is defined the same way from $(x^i_f)$.

Formalization targets

Goal: Theorem 3.2 (p. 840)

Let $z^*$ be a fixed point of $T$, $x^* = \operatorname{prox}_{\gamma g}(z^*)$, and suppose $f$ is $L$-Lipschitz continuous on the closed ball $B\big(x^*, (1+\gamma/\beta)\|z^0 - z^*\|\big)$. Then there is a constant $C$, independent of $k$, with

$$(f+g+h)(\bar x^k_g) - (f+g+h)(x^*) \le \frac{C}{k+1}\qquad (k \ge 0).$$

The goal asserts the order $O(1/(k+1))$ and leaves the constant free, so it is not invalidated by a sharper constant.

Milestones

1. Corollary 2.1, Part 1 (p. 834): $\|z^j - z^*\|$ is nonincreasing.
2. Lemma 3.1 (p. 838): $x^j_f, x^j_g \in B\big(x^*, (1+\gamma/\beta)\|z^0 - z^*\|\big)$ for all $j$.
3. Eq. (3.2) (p. 839): for all $k \ge 0$,

$$2\gamma\big(f(x^k_f) + g(x^k_g) + h(x^k_g) - (f+g+h)(x^*)\big) \le \|z^k - x^*\|^2 - \|z^{k+1} - x^*\|^2 - \|z^k - z^{k+1}\|^2 + 2\gamma\langle z^k - z^{k+1}, \nabla h(x^k_g)\rangle .$$

4. Theorem 3.1 (p. 838): the last-iterate rate $(f+g+h)(x^k_g) - (f+g+h)(x^*) = o\big(1/\sqrt{k+1}\big)$.
5. Eq. (2.7) (p. 836), with $\lambda_k \equiv 1$: for $\gamma/(2\beta) < \varepsilon < 1$,

$$\sum_{i=k}^\infty \|\nabla h(x^i_g) - \nabla h(x^*)\|^2 \le \frac{\|z^k - z^*\|^2}{\gamma(2\beta - \gamma/\varepsilon)} .$$

6. Eq. (3.4) (p. 840): $\|\bar x^k_f - \bar x^k_g\| \le 5\|z^0 - z^*\|/(k+1)$.

Significance

The result. Theorem 3.1 gives the last iterate an objective error of $o(1/\sqrt{k+1})$. Theorem 3.2 shows that averaging with linearly increasing weights improves this to $O(1/(k+1))$, the rate of the standard uniform ergodic average, while putting more weight on recent iterates. The paper notes that this matters when the iterates $x^k_g$ are sparse vectors or low-rank matrices and the average should stay close to them. The rates hold under a local Lipschitz condition on one of the two nonsmooth terms only, so $g$ may be the indicator function of a constraint set. They therefore cover the constrained applications of Section 4 of the paper.

Formalizing it. The results are proved in the paper. No machine-checked version of this scheme or its rates exists on the platform or, as far as is known, in Mathlib. A formalization produces a checked proof in an arbitrary real Hilbert space with extended-valued $f, g$. It also produces infrastructure that Mathlib lacks: proximal maps characterised by minimisation, the prox-subgradient inclusion, Fejér monotonicity of an averaged-operator iteration, and a weighted Jensen inequality for extended-valued convex functions. All of these can be reused by other splitting and proximal-gradient missions. The formalization also checks the constants: the last display of the published proof of Theorem 3.2 drops a factor $2\gamma$ in front of the Lipschitz term, and the printed ball in both theorems is centred at $0$ where the proof needs $x^*$.

Difficulty

The obvious argument sums the one-step inequality (3.2). That controls the objective at the two different points $x^k_f$ and $x^k_g$, and only $f(x^k_f)$ appears, never $f(x^k_g)$. Moving from one point to the other needs the Lipschitz hypothesis on $f$, and so it needs every iterate, and every weighted average, to stay in the ball on which that hypothesis holds. For the weighted average there is a further obstacle: the cross term $2\gamma\langle z^k - z^{k+1}, \nabla h(x^k_g)\rangle$ does not telescope under the weights $(i+1)$. Controlling it requires the summability of the gradient differences (2.7), which is inherited from the averagedness analysis of Section 2 and not from convexity alone. Uniform averaging with the same argument does not give the weighted statement, and the weights must not be replaced.

Formalization scope

• $H$ is an arbitrary real Hilbert space (InnerProductSpace ℝ H, CompleteSpace H), not $\mathbb R^n$.
• $f, g : H \to$ EReal. They are proper (never $\bot$, somewhere $\ne \top$), lower semicontinuous, and have a convex epigraph in $H \times \mathbb R$. $h : H \to \mathbb R$ is convex and differentiable, and Mathlib's gradient h is $\beta^{-1}$-Lipschitz.
• Proximal maps are not constructed. A map $P$ is assumed to minimise $f(y) + \|y - x\|^2/(2\gamma)$ for every $x$. Such a map exists and is unique for closed proper convex $f$, so nothing is lost.
• Algorithm 2 is fixed with $\lambda_k \equiv 1$, the only case of Theorems 3.1 and 3.2. Iterates are indexed from $0$. The fixed point $z^*$ is a hypothesis, $T z^* = z^*$, and $x^* := \operatorname{prox}_{\gamma g}(z^*)$. Assumption 1 of the paper follows from this and is not assumed separately.
• Ball centre. The theorems print $B(0, (1+\gamma/\beta)\|z^0 - z^*\|)$. The proofs use Lemma 3.1, whose ball is centred at $x^*$, so the ball here is centred at $x^*$. "$f$ is $L$-Lipschitz on the ball" is stated as: $f$ is finite on the ball, and its real-valued restriction is $L$-Lipschitz there.
• O(·) and o(·). $O(1/(k+1))$ is $\exists C \in \mathbb R,\ \forall k,\ (f+g+h)(\bar x^k_g) \le (f+g+h)(x^*) + C/(k+1)$, with $C$ chosen after all data. $o(1/\sqrt{k+1})$ is $\sqrt{k+1}\,\big((f+g+h)(x^k_g) - (f+g+h)(x^*)\big) \to 0$, together with finiteness of the objective values as part of the conclusion. No explicit constant from the proof is stated, because the published constant drops a factor.
• Corollary 2.1 Part 1 and Eq. (2.7) are stated for Algorithm 2 with $\lambda_k \equiv 1$, $\gamma \in (0, 2\beta)$ and $\varepsilon \in (\gamma/(2\beta), 1)$. As printed, Corollary 2.1's condition on $\tau_k$ excludes $\lambda_k \equiv 1$, but Section 3 uses Part 1 in exactly this case. Summability in (2.7) is part of the conclusion.
• Trivialization ruled out. Objective values are extended reals, and the goal compares them without subtraction. The value $(f+g+h)(x^*)$ is proved finite as part of the conclusion. So the goal cannot hold through $\infty - \infty$ or through an infinite right-hand side.

Welcome contributions: the prox–subgradient inclusion for EReal-valued convex functions, averagedness and Fejér monotonicity of $T$ (the companion mission on Section 2 treats the general operator case), a weighted Jensen inequality in EReal, and proofs of the milestones in the listed order.

Selected references

• D. Davis and W. Yin, A Three-Operator Splitting Scheme and its Optimization Applications, Set-Valued and Variational Analysis 25 (2017) 829–858. https://doi.org/10.1007/s11228-017-0421-z
• H. H. Bauschke and P. L. Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Spaces, 2nd ed., Springer, 2017. https://doi.org/10.1007/978-3-319-48311-5
• 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
• D. Davis and W. Yin, Convergence Rate Analysis of Several Splitting Schemes, in Splitting Methods in Communication, Imaging, Science, and Engineering, Springer, 2016. https://doi.org/10.1007/978-3-319-41589-5_4

Motivation

Many problems in convex optimization, variational inequalities and signal processing reduce to a monotone inclusion: find a point $x$ at which the sum of several monotone operators contains $0$. When the sum has two terms, the classical operator-splitting methods (Douglas–Rachford, forward–backward, forward–backward–forward) solve it by iterating a fixed-point map that uses each operator separately, through its resolvent or through a forward (explicit) step. Problems with three terms, for instance a smooth loss plus two nonsmooth regularizers or constraints, are common in practice, and before 2015 no fixed-point map was known that handled three operators one at a time without a product-space reformulation.

Davis and Yin (Set-Valued Var. Anal. 25 (2017) 829–858; preprint arXiv:1504.01032) introduced such a map, now called Davis–Yin three-operator splitting. It contains Douglas–Rachford splitting ($C = 0$) and forward–backward splitting ($B = 0$) as special cases, and it has become a standard building block of first-order methods for composite optimization. This mission formalizes Section 2 of the paper: the fixed-point encoding, the averagedness of the map, and the weak and strong convergence of the resulting iteration.

Setting

Let $H$ be a real Hilbert space. A set-valued operator $A : H \to 2^H$ is monotone if $\langle x - y, u - v\rangle \ge 0$ for all $u \in Ax$, $v \in Ay$, and maximal monotone if its graph is not properly contained in the graph of another monotone operator. Its domain is $\operatorname{dom}(A) = \{x : Ax \ne \emptyset\}$ and the zero set of an operator $M$ is $\operatorname{zer}(M) = \{x : 0 \in Mx\}$. A single-valued $C : H \to H$ is $\beta$-cocoercive ($\beta > 0$) if $\beta\|Cx - Cy\|^2 \le \langle Cx - Cy, x - y\rangle$ for all $x, y$.

Problem (1.1) is: given maximal monotone $A, B$ and $\beta$-cocoercive $C$, find

$$x \in H \quad\text{with}\quad 0 \in Ax + Bx + Cx .$$

For $\gamma > 0$ the resolvent $J_{\gamma A} = (I + \gamma A)^{-1}$ is the map with $x \in J_{\gamma A}x + \gamma A(J_{\gamma A}x)$. The Davis–Yin operator (Eq. (1.2)) is

$$T := J_{\gamma A} \circ (2J_{\gamma B} - I - \gamma C \circ J_{\gamma B}) + I - J_{\gamma B}.$$

Algorithm 1 starts from $z^0 \in H$ and, for relaxation parameters $\lambda_k > 0$, iterates

$$x_B^k = J_{\gamma B}(z^k),\qquad x_A^k = J_{\gamma A}(2x_B^k - z^k - \gamma Cx_B^k),\qquad z^{k+1} = z^k + \lambda_k(x_A^k - x_B^k),$$

so that $z^{k+1} = (1 - \lambda_k)z^k + \lambda_k Tz^k$. A sequence converges weakly, $u_k \rightharpoonup u$, if $\langle u_k, y\rangle \to \langle u, y\rangle$ for every $y \in H$.

Formalization targets

Goal: Theorem 2.1 (Main convergence theorem)

Fix $\varepsilon \in (0,1)$, $\gamma \in (0, 2\beta\varepsilon)$, $\alpha = 1/(2-\varepsilon)$ and $\lambda_k \in (0, 1/\alpha)$ with $\sum_k \tau_k = \infty$, where $\tau_k = \lambda_k(1-\lambda_k) + \lambda_k(1-\alpha)/\alpha$, and $\inf_k \lambda_k > 0$. If $\operatorname{Fix} T \ne \emptyset$, there is $z^* \in \operatorname{Fix} T$ with $z^k \rightharpoonup z^*$ and

$$Cx_B^k \to Cx^* \ \ (\forall x^* \in \operatorname{zer}(A+B+C)),\qquad x_B^k \rightharpoonup J_{\gamma B}(z^*) \in \operatorname{zer}(A+B+C),\qquad x_A^k \rightharpoonup J_{\gamma B}(z^*),$$

and if $A$ or $B$ is uniformly monotone on every nonempty bounded subset of its domain, or $C$ is demiregular at every zero of $A + B + C$, then $x_B^k$ and $x_A^k$ converge strongly to a common point of $\operatorname{zer}(A + B + C)$.

Milestones

In the order the proof uses them: Lemma 2.1 (the identities for one application of $T$), Lemma 2.2 ($\operatorname{zer}(A+B+C) = J_{\gamma B}(\operatorname{Fix} T)$), Lemma 2.3 (inequality (2.1)), Proposition 2.1 ($T$ is $2\beta/(4\beta-\gamma)$-averaged, inequality (2.2)), Remark 2.1 (the strengthened inequality (2.4)), Corollary 2.1 Parts 1–3 (Fejér monotonicity, vanishing residual, weak convergence of $z^k$), Corollary 2.1 Part 4 (the residual rates $\|Tz^k - z^k\|^2 \le \|z^0 - z^*\|^2/(\underline\tau(k+1))$ and $o(1/(k+1))$), and Eqs. (2.6)–(2.7) (the per-step descent inequality and its summed form).

Significance

Theorem 2.1 is the basic convergence guarantee for three-operator splitting: it certifies that the computable sequences $x_B^k$, $x_A^k$, not only the auxiliary sequence $z^k$, approach a solution of (1.1). In infinite dimensions this is the delicate part: for Douglas–Rachford splitting ($C = 0$) weak convergence of the shadow sequence $J_{\gamma B}(z^k)$ was only established by Svaiter in 2011. The result underlies the convergence of the many algorithms obtained from it by specialization (Douglas–Rachford, forward–backward, and the three-block methods of Section 4 of the paper), and the averagedness coefficient of Proposition 2.1 reduces, for $B = 0$, to the best known one for forward–backward splitting.

All statements of this mission are proved in the paper, partly by appeal to Bauschke and Combettes' monograph (Krasnosel'skiĭ–Mann convergence, the demiclosedness of maximal monotone graphs). None of them has a machine-checked proof: Mathlib has no maximal monotone operators, resolvents, averaged maps or Krasnosel'skiĭ–Mann theorem. The mission therefore produces both a formal proof of the Davis–Yin theorem and a first body of monotone-operator theory in Lean.

Difficulty

The fixed-point part is standard once $T$ is known to be averaged: Krasnosel'skiĭ–Mann theory and Opial's argument give $z^k \rightharpoonup z^*$. The obstacle is transferring this to $x_B^k = J_{\gamma B}(z^k)$. Resolvents are nonexpansive but not weakly continuous, so $z^k \rightharpoonup z^*$ does not imply $J_{\gamma B}(z^k) \rightharpoonup J_{\gamma B}(z^*)$; the naive argument fails at exactly this step. Identifying the weak cluster points of $x_B^k$ requires a closedness property of sums of maximal monotone operators under mixed weak and strong convergence, fed by the strong convergence of $Cx_B^k$, which in turn needs the extra term of (2.4) that (2.2) discards. Strong convergence in Part 2 needs yet another argument for each of the three alternative hypotheses.

Formalization scope

• $H$ is an arbitrary real Hilbert space (NormedAddCommGroup, InnerProductSpace ℝ, CompleteSpace); a finite-dimensional space would identify weak and strong convergence and change the theorems.
• Operators $A, B$ are H → Set H; $C$ is single-valued H → H. The resolvents are not constructed: $J_A, J_B$ are maps satisfying the resolvent inclusion $\gamma^{-1}(x - Jx) \in A(Jx)$, which for maximal monotone operators determines them uniquely and exists by Minty's theorem.
• Weak convergence is $\langle u_k, y\rangle \to \langle u, y\rangle$ for every $y$; strong convergence is norm convergence. Iterates are indexed from $0$.
• The printed hypothesis $\alpha = 1/(2-\varepsilon) < 2\beta/(4\beta-\gamma)$ of Corollary 2.1 and Theorem 2.1 contradicts $\gamma < 2\beta\varepsilon$ (it is a typo for $>$) and is not assumed. The printed $\tau_k = (1-\lambda_k/\alpha)\lambda_k/\alpha$ is replaced by the $\tau_k$ of the proof (p. 836), a weaker hypothesis.
• Uniform monotonicity uses a nondecreasing $\varphi : [0,\infty) \to [0,+\infty]$ with $\varphi(0) = 0$ that vanishes only at $0$, as the proof requires; with $\varphi \equiv 0$ allowed, Part 2(a) would be false.
• The O-constant of Corollary 2.1 Part 4 is explicit, $\|z^0 - z^*\|^2/\underline\tau$, and the little-$o$ is stated as $(k+1)\|Tz^k - z^k\|^2 \to 0$. Eq. (2.7) is stated with a uniform lower bound $\underline\lambda \le \lambda_i$ in place of the printed $\lambda_k$, with summability part of the conclusion.
• A formalization with $T$ an arbitrary averaged map, with resolvents replaced by arbitrary nonexpansive maps, or with the contradictory comparison of $\alpha$ kept as a hypothesis would make the theorem vacuous or different; all three are ruled out.

A complete development needs the basic theory of monotone operators (monotonicity of resolvents' graphs, firm nonexpansiveness of resolvents, weak-to-strong closedness of maximal monotone graphs), Krasnosel'skiĭ–Mann iteration with Opial's lemma, and weak sequential compactness of bounded sets in Hilbert space. All of this is reusable far beyond this mission, and contributions of any of these pieces as separate theorems are welcome.

Selected references

• D. Davis and W. Yin, A Three-Operator Splitting Scheme and its Optimization Applications, Set-Valued and Variational Analysis 25 (2017) 829–858. https://doi.org/10.1007/s11228-017-0421-z
• H. H. Bauschke and P. L. Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Spaces, Springer, 2011. https://doi.org/10.1007/978-1-4419-9467-7
• 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
• D. Davis and W. Yin, Convergence rate analysis of several splitting schemes, in: Splitting Methods in Communication, Imaging, Science, and Engineering, Springer, 2016. https://arxiv.org/abs/1406.4834

Motivation

In 1955 E. M. L. Beale showed how Dantzig's simplex method, which was built for linear objectives, can be carried over to certain nonlinear convex objectives that are minimized subject to linear inequalities (Beale 1955). Section 4 of that paper treats one such objective: the sum of the $t$ largest of a set of $g$ linear forms. Beale's motivation comes from the theory of games: "if the enemy has to choose $t$ out of a set of $g$ possible actions, and $L_f$ represents his average gain through using the $f$th", then the defender wants to minimize the sum of the $t$ largest $L_f$.

The same objective can be written as a linear program. One introduces a bound $u$ and requires every sum of $t$ forms to be at most $u$. That formulation has $\binom{g}{t}$ constraints, which is unwieldy once $t>1$ and $g$ is large. Beale's alternative works with the nonlinear objective directly, and he needs a test that tells him when the current basic solution is already optimal. This mission formalizes that test, Theorem 1 of the paper.

The objective reappears in later work under other names: the sum of the $k$ largest components of a vector, the "top-$k$ sum", and $k$ times the conditional value-at-risk of an empirical distribution. Beale's paper is an early source for its optimality conditions.

Setting

There are real variables $z_l$, indexed by $l$ in a finite set (possibly empty), and $u_1,\dots,u_s$. Two linear forms in these variables are given,

$$A=A_0+\sum_l A_l z_l+\sum_{f=1}^{s}\varphi_f u_f,\qquad L_0=c_{00}+\sum_l c_{0l} z_l+\sum_{f=1}^{s}\theta_f u_f,$$

together with $s$ further forms

$$L_f=L_0-u_f\qquad(f=1,\dots,s).$$

For an integer $\tau\ge0$ the objective is

$$C=A+\bigl(\text{sum of the }\tau\text{ largest of }L_0,L_1,\dots,L_s\bigr).$$

The sum of the $\tau$ largest of $s+1$ numbers is the largest total of any $\tau$ of them. Ties do not make it ambiguous.

The feasible region is fixed by a set $F$ of indices. The variables $z_l$ with $l\in F$ and all the $u_f$ are free, and every other $z_l$ is restricted to $z_l\ge0$. At the origin $z=0$, $u=0$ all $s+1$ forms are equal to $c_{00}$, so the origin is where $C$ fails to be differentiable. In Beale's algorithm the origin is the current basic solution: the $u_f$ measure how far the "borderline" forms sit from a chosen critical form, and $A$ collects the forms that are certainly among the largest.

Write $a_l=A_l+\tau c_{0l}$ and $w_f=\varphi_f+\tau\theta_f$.

Formalization targets

Goal: Theorem 1 (a), p. 179

For $\tau\le s$, $C$ is minimized over the feasible region when all the $z_l$ and $u_f$ vanish if and only if

$$\begin{aligned} &a_l\ge0\ \text{for all } l, \qquad a_l=0\ \text{for all } l\in F,\\ &0\le w_f\le1\ \text{for all } f,\qquad \tau-1\le\sum_{f=1}^{s}w_f\le\tau . \end{aligned}\tag{4.5}$$

"Minimized" means a global minimum: $C(0,0)\le C(z,u)$ at every feasible point.

Milestones

1. Convexity (p. 179). $C$ is a convex function of $(z,u)$ for $\tau\le s+1$.
2. Descent rules (second half of Theorem 1 (a), p. 179). When a condition of (4.5) fails, a stated move of one variable, or of all $u_f$ together, lowers $C$ below $C(0,0)$ for every small enough step. There are six moves: $z_l\uparrow$ if $a_l<0$; $z_l\downarrow$ if $a_l>0$ and $l\in F$; $u_f\uparrow$ if $w_f<0$; $u_f\downarrow$ if $w_f>1$; all $u_f\uparrow$ if $\sum w_f<\tau-1$; all $u_f\downarrow$ if $\sum w_f>\tau$.
3. The rearrangement identity (proof of Theorem 1 (a), p. 180). If $1\le\tau\le s$, $u'_1\le\dots\le u'_s$ and $u'_\tau\le0$, then

$$C=A_0+\tau c_{00}+\sum_l a_l z'_l+\sum_{f=1}^{\tau}(w_f-1)(u'_f-u'_\tau)+\sum_{f=\tau+1}^{s}w_f(u'_f-u'_\tau)+\Bigl\{\sum_{f=1}^{s}w_f-\tau\Bigr\}u'_\tau .$$

4. Theorem 1 (b) (p. 180). For $\tau=s+1$, the origin is a minimum if and only if (4.5) holds and $w_f=1$ for every $f$. Otherwise some value of $u_f$ with the sign opposite to $w_f-1$ lowers $C$.

Significance

Theorem 1 is the optimality test of Beale's simplex method for the sum-of-largest objective. The algorithm on pp. 178–179 changes nonbasic variables one at a time. When no single change is profitable it applies Theorem 1: either (4.5) holds and the current solution is optimal, or one of the six descent rules names the variable to change next. The test is exact even though the objective is not differentiable at the current point. It is a closed-form description of the subdifferential of a top-$\tau$ sum at a point where all the forms tie. The theorem is also the base case of the multi-group generalization that Beale mentions on p. 181.

The paper proves Theorem 1 by hand. To our knowledge neither the theorem nor the rearrangement identity behind it has been formalized in any proof assistant. The mission produces:

• a checked statement and proof of the test, including the degenerate cases $\tau=0$ and $s=0$, which the paper does not discuss separately;
• the boundary case $\tau=s+1$;
• a reusable Lean definition of the sum of the $\tau$ largest entries of a finite real family, with its convexity.

Difficulty

Necessity, the "only if" direction, is the part the paper calls obvious: each descent rule changes $C$ linearly for small steps. Two features still have to be handled explicitly. The step must be small only in rule-dependent ways, and the ordering of the forms changes along the moves of rules 4 and 6.

Sufficiency is where the work lies. The naive argument, "the directional derivative in every coordinate direction is non-negative, so the origin is a minimum", fails because $C$ is not differentiable at the origin. Nonnegative derivatives along the coordinate axes do not control mixed directions in which several $u_f$ move by different amounts, which reorders the forms. Which $\tau$ forms are the largest then depends on the point, and the paper settles the configurations in which $L_0$ is among the $\tau$ largest by an informal appeal to the "essential symmetry" between $L_0$ and the other forms. A formal proof cannot leave that appeal informal: the forms are parametrised relative to $L_0$ (each $L_f$ is $L_0-u_f$), so the symmetry is a change of variables that has to be written down and shown to preserve (4.5).

Formalization scope

• Data. The variables are z : Fin r → ℝ (any r, including $0$) and u : Fin s → ℝ. The paper's $u_f$ for $f=1,\dots,s$ is Lean's u f for $f=0,\dots,s-1$. The coefficients $(A_0,A_l,\varphi_f,c_{00},c_{0l},\theta_f)$ form a structure Forms r s.
• Forms. The family $L_0,\dots,L_s$ is Fin (s+1) → ℝ, with index $0$ for $L_0$ and index f.succ for $L_0-u_f$. The free set $F$ is a Finset (Fin r), and $\tau$ is a natural number cast to $\mathbb R$ wherever it multiplies a coefficient.
• Sum of the largest. sumLargest τ v is the maximum over $\tau$-element subsets $S$ of $\sum_{i\in S}v_i$ (Finset.sup' over powersetCard). It is the junk $0$ for $\tau$ larger than the number of entries, a case no statement uses.
• Minimality. "Minimized when all variables vanish" is the global statement $C(0,0)\le C(z,u)$ for all $(z,u)$ with $z_l\ge0$ for $l\notin F$. It is not a local minimum, and the sign constraints on restricted $z_l$ are kept: they are why the first condition of (4.5) is an inequality.
• Descent. "$C$ can be decreased by moving $x$ from zero" is a strict decrease for all step sizes in some interval $(0,\varepsilon)$, with every other variable at zero.
• No trivialization. The goal is an equivalence with no hypothesis beyond $\tau\le s$. Neither direction can be satisfied vacuously, and the cases $\tau=0$ and $s=0$ are included, as on the page.
• Added hypotheses. The rearrangement milestone assumes $\tau\ge1$, because the paper's $u'_\tau$ does not exist at $\tau=0$. Its second line uses $c_{0l}$ where the page misprints $c_l$.

Needed infrastructure:

• basic lemmas on sumLargest: its value at a constant family, at a family sorted by a monotone shift, and under adding a common constant;
• the change of variables behind the paper's symmetry between $L_0$ and the other forms.

These lemmas are reusable for any top-$k$-sum or empirical-CVaR objective. Contributions are welcome at any level: lemmas about sumLargest, any of the milestones, or an alternative sufficiency proof through convexity and one-sided directional derivatives.

Not in scope: the pivoting rules (4.2)–(4.4), the degeneracy discussion on pp. 180–181, and the multi-group generalization, which the paper says is "cumbersome to state" and does not state.

Selected references

• E. M. L. Beale, On Minimizing a Convex Function Subject to Linear Inequalities, Journal of the Royal Statistical Society, Series B 17(2), 173–184, 1955. https://doi.org/10.1111/j.2517-6161.1955.tb00191.x
• G. B. Dantzig, A. Orden and P. Wolfe, The generalized simplex method for minimizing a linear form under linear inequality restraints, Pacific Journal of Mathematics 5(2), 183–195, 1955. https://doi.org/10.2140/pjm.1955.5.183
• R. T. Rockafellar and S. Uryasev, Optimization of conditional value-at-risk, Journal of Risk 2(3), 21–41, 2000. https://doi.org/10.21314/JOR.2000.038

Motivation

Interior point methods for a linear program $\min\{c^\top x : Ax\ge b\}$ with $A\in\mathbb R^{m\times n}$ follow the central path of the logarithmic barrier $-\sum_i\log s_i$, where $s=Ax-b$ is the slack vector. Renegar's path-following analysis (1988) gives $O(\sqrt m\,L)$ iterations, and for decades this was the best bound for methods whose iterations cost a linear system solve. Vaidya's volumetric barrier $-\log\det(A^\top S^{-2}A)$ and the hybrid volumetric barriers of Vaidya and of Anstreicher (references [45] and [2] of the paper) reached $O((m\,\mathrm{rank}(A))^{1/4}L)$ iterations at the price of more expensive linear algebra. Nesterov and Nemirovski showed that a universal barrier gives $O(\sqrt n\,L)$ iterations, but that barrier cannot be evaluated efficiently.

Lee and Sidford (FOCS 2014; full version arXiv:1312.6677) obtained $\tilde O(\sqrt{\mathrm{rank}(A)}\,L)$ iterations, each costing $\tilde O(1)$ linear system solves, by following a weighted central path whose weights are recomputed from the slacks. The weights come from a weight function $g$, defined as the minimizer of a regularized D-optimal-design problem. This mission is about that weight function and the theorem (Theorem 1 of the paper) certifying its properties. The companion mission, Path-Finding Methods for Linear Programming I, formalizes the path-following framework (Theorem 5 of §IV.C) that consumes these properties.

Setting

Fix $A\in\mathbb R^{m\times n}$ with full column rank, $\mathrm{rank}(A)=n$, and $1\le n<m$. For vectors $s,w\in\mathbb R^m_{>0}$ write $S=\mathrm{diag}(s)$, $W=\mathrm{diag}(w)$, $W^\alpha=\mathrm{diag}(w_i^\alpha)$, and $A_s=S^{-1}A$. For a matrix $M$ let $\|v\|_M=\sqrt{v^\top Mv}$.

Projection matrix and slack sensitivity (Definition 2, p. 428). The projection matrix is $P_{S^{-1}A}(w)=W^{1/2}S^{-1}A\,(A^\top S^{-1}WS^{-1}A)^{-1}A^\top S^{-1}W^{1/2}$, and the slack sensitivity is

$$\gamma(s,w)=\max_{i\in[m]}\big\|W^{-1/2}\mathbb 1_i\big\|_{P_{S^{-1}A}(w)} .$$

Weight function (Definition 4, p. 428). A map $g:\mathbb R^m_{>0}\to\mathbb R^m_{>0}$ is a weight function with constants $c_1,c_\gamma,c_r$ if it is differentiable and, for every $s>0$, with $G(s)=\mathrm{diag}(g(s))$, $G'(s)$ the Jacobian of $g$ at $s$, and $\|y\|_{G(s)}=\sqrt{\sum_ig_i(s)y_i^2}$:

1. Size: $\|g(s)\|_1\le c_1$;
2. Slack sensitivity: $c_\gamma\ge1$ and $\gamma(s,g(s))\le c_\gamma$;
3. Step consistency: $c_r\ge1$ and for all $r\ge c_r$, $y\in\mathbb R^m$: $\|(I+r^{-1}G^{-1}G'S)y\|_{G(s)}\le\|y\|_{G(s)}$ and $\|y+r^{-1}G^{-1}G'Sy\|_\infty\le\|y\|_\infty+c_r\|y\|_{G(s)}$;
4. Uniformity: $\|g(s)\|_\infty\le2$.

The regularized objective (6), p. 429. For $\alpha,\beta\in\mathbb R$,

$$\hat f(s,w)=\mathbb 1^\top w-\frac1\alpha\log\det\big(A_s^\top W^\alpha A_s\big)-\beta\sum_{i\in[m]}\log w_i ,\qquad g(s)=\arg\min_{w\in\mathbb R^m_{>0}}\hat f(s,w).$$

At $\alpha=1,\beta=0$ this is the D-optimal design problem, dual to computing the John ellipsoid of the polytope $\{y:|[A(y-x)]_i|\le s_i\}$ (§V.B).

Formalization targets

Goal: Theorem 1 (Properties of Weight Function), §V.A, p. 429

With

$$\alpha=1-\Big(\log_2\frac{2m}{\mathrm{rank}(A)}\Big)^{-1},\qquad \beta=\frac{\mathrm{rank}(A)}{2m},$$

the objective $\hat f(s,\cdot)$ has a unique minimizer over $\mathbb R^m_{>0}$ for every $s>0$, and the resulting $g$ is a weight function with

$$c_1(g)=2\,\mathrm{rank}(A),\qquad c_\gamma(g)=2,\qquad c_r(g)=2\log_2\frac{2m}{\mathrm{rank}(A)} .$$

Milestones: the three bullets of Theorem 1

• Size: every minimizer $w$ of $\hat f(s,\cdot)$ satisfies $\|w\|_1\le2\,\mathrm{rank}(A)$.
• Slack sensitivity: every minimizer $w$ satisfies $\gamma(s,w)\le2$.
• Step consistency: any map $g$ selecting a minimizer at every $s>0$ is differentiable on $\mathbb R^m_{>0}$ and satisfies the two step-consistency inequalities for every $r\ge2\log_2\frac{2m}{\mathrm{rank}(A)}$.

A supporting (non-milestone) item states the existence and uniqueness of the minimizer on its own.

Significance

The result. Theorem 1 is the input that turns the weighted path-following framework into an $\tilde O(\sqrt{\mathrm{rank}(A)}\,L)$-iteration method: the framework needs $O(c_\gamma^{-1}c_r^{-3}c_1^{-1/2})$-sized steps in $t$ (p. 428), and Theorem 1 makes that $\tilde\Omega(1/\sqrt{\mathrm{rank}(A)})$. The step consistency bound is what allows the weights to be recomputed after each Newton step without losing centrality. The same construction underlies later work on Lewis-weight barriers and on fast approximate John ellipsoids and maximum flow (§VIII of the paper).

Formalizing it. The theorem is proved in the full version of the paper (arXiv:1312.6677); the FOCS extended abstract contains no proofs. No part of it has a machine-checked proof. A complete formalization would give a verified account of leverage-score calculus (sums of leverage scores equal the rank; derivatives of projection matrices), of the convexity of $w\mapsto-\log\det(A^\top W^\alpha A)$ for $\alpha\in(0,1)$, and of differentiability of an argmin via the implicit function theorem, none of which is currently packaged in Mathlib in this form.

Difficulty

Size and slack sensitivity are statements about the minimizer, which is only characterized implicitly; they require precise matrix calculus for $\log\det(A_s^\top W^\alpha A_s)$ and a comparison between the matrices $A^\top WA$ (which defines $\gamma$) and $A^\top W^\alpha A$ (which defines $g$). The specific values of $\alpha$ and $\beta$ matter here: the unregularized choice $\alpha=1$, $\beta=0$ makes the problem degenerate (p. 429).

The hard part is step consistency. The Jacobian $G'$ of an argmin is available only implicitly, as the solution of a linear system obtained by differentiating the optimality condition. A bound on $\|G'\|$ that depends on $m$ is easy to get and useless: the theorem needs the operator norm of $I+r^{-1}G^{-1}G'S$ in the $G(s)$-norm to be at most $1$ as soon as $r$ exceeds $2\log_2(2m/\mathrm{rank}(A))$, and an $\ell_\infty$ bound with only an additive $c_r\|y\|_{G(s)}$ loss.

Existence and differentiability of the minimizer are conclusions, not hypotheses. The minimization is over an open orthant on which the objective is not obviously coercive or strictly convex for $\alpha<1$, and differentiability of $g$ requires the Hessian of $\hat f$ at the minimizer to be invertible.

Formalization scope

Vectors are Fin m → ℝ, matrices Matrix (Fin m) (Fin n) ℝ; inverses are Matrix.inv, $\log\det$ is Real.log (Matrix.det …), $w_i^\alpha$ is Real.rpow, $\log_2$ is Real.logb 2, the Jacobian is fderiv ℝ g s, and $\|\cdot\|_\infty$ is Mathlib's sup norm on Fin m → ℝ.

Conventions and pinned hypotheses:

• Full column rank A.rank = n is assumed in every theorem. The paper never states it, but without it $A_s^\top W^\alpha A_s$ is singular and every formula is undefined (in Lean, Matrix.inv and Real.log would return junk $0$).
• $1\le n<m$. $\beta=\mathrm{rank}(A)/(2m)$ and $\log_2(2m/\mathrm{rank}(A))$ need $\mathrm{rank}(A)\ge1$; at $m=\mathrm{rank}(A)$ the page's $\alpha$ is $0$ and $1/\alpha$ in (6) is undefined.
• Reading of $\alpha$: the exponent $-1$ is the reciprocal of $\log_2\frac{2m}{\mathrm{rank}(A)}$, giving $\alpha\in(0,1)$.
• Size is an upper bound $\|g(s)\|_1\le c_1$ (the paper's weight function has $\|g(s)\|_1=\tfrac32\mathrm{rank}(A)$, while Theorem 1 reports $c_1=2\,\mathrm{rank}(A)$).
• The first step-consistency bullet (an operator-norm bound) is stated for every vector $y$.
• $g$ is any map $\mathbb R^m\to\mathbb R^m$ whose value at each positive $s$ minimizes $\hat f(s,\cdot)$ over $\mathbb R^m_{>0}$. Only its values on the open orthant matter. The goal also asserts that such minimizers exist and are unique, so it is not vacuous.

Ruling out trivializations: the goal does not assume $g$ to be a weight function or to be differentiable, and it does not replace $g$ by an arbitrary weight function; differentiability is a conclusion (a predicate using fderiv without it would make step consistency hold vacuously wherever $g$ fails to be differentiable).

Useful infrastructure, reusable beyond this mission: leverage scores and their sum; derivatives of $w\mapsto\log\det(A^\top WA)$ and of projection matrices; convexity of $-\log\det(A^\top W^\alpha A)$ in $w$ (related to the published ConvexOptimization.log_det_concaveOn); differentiability of the argmin of a strictly convex smooth function. Contributions of these as separate theorems are welcome, as is a proof of any single bullet of Theorem 1.

Selected references

• Y. T. Lee, A. Sidford, Path Finding Methods for Linear Programming: Solving Linear Programs in Õ(√rank) Iterations and Faster Algorithms for Maximum Flow, FOCS 2014, pp. 424–433. https://doi.org/10.1109/FOCS.2014.52
• Y. T. Lee, A. Sidford, Path Finding I: Solving Linear Programs with Õ(√rank) Linear System Solves, arXiv, 2013. https://arxiv.org/abs/1312.6677
• J. Renegar, A polynomial-time algorithm, based on Newton's method, for linear programming, Mathematical Programming 40 (1988). https://doi.org/10.1007/BF01580724