Numerical Analysis

Motivation

Obstacle problems, contact problems and control-constrained optimal control problems lead, after discretization, to complementarity systems of the form "find $y \le \psi$ and $\lambda \ge 0$ with $Ay + \lambda = f$ and $\lambda_i(y_i - \psi_i) = 0$". The primal-dual active set strategy solves such systems by guessing which constraints are active, solving one linear system, and updating the guess. It is widely used because each iteration is a single linear solve and, in practice, it terminates after few iterations.

Hintermüller, Ito and Kunisch (SIAM J. Optim. 13(3), 2002; authors' version hal-01660511v1) showed that the strategy is a semismooth Newton method applied to a reformulation of the complementarity system with a max-function. This gives local superlinear convergence. For global convergence, from arbitrary initial data, the paper gives three sufficient conditions on the matrix $A$: the M-matrix case (Theorem 3.2), the condition of Theorem 3.3 treated in this mission, and small perturbations of M-matrices (Theorem 3.4). Theorem 3.3 is the one that goes beyond sign patterns: it covers P-matrices whose blocks satisfy a norm bound and a cone condition, with the sum of the components of $y$ serving as a merit function.

Setting

Fix $n \ge 0$, a matrix $A \in \mathbb{R}^{n \times n}$, vectors $f, \psi \in \mathbb{R}^n$ and a parameter $c > 0$. Vectors are ordered componentwise. The problem is

$$Ay + \lambda = f, \qquad \lambda - \max\big(0, \lambda + c(y - \psi)\big) = 0, \tag{3.1}$$

with the maximum taken componentwise. In the paper's notation the matrix is $A$ and the index sets are calligraphic: $\mathcal{A}$ (active) and $\mathcal{I}$ (inactive).

A P-matrix is a square matrix all of whose principal minors are positive. For an index set $\mathcal{I} \subseteq \{1,\dots,n\}$ with complement $\mathcal{A}$, $A_{\mathcal{I}}$ is the principal submatrix of $A$ on $\mathcal{I}$, and $A_{\mathcal{I}\mathcal{A}}$ is the submatrix with rows in $\mathcal{I}$ and columns in $\mathcal{A}$. For a P-matrix every $A_{\mathcal{I}}$ is invertible.

The primal-dual active set algorithm starts from an arbitrary pair $(y^0, \lambda^0)$. Given $(y^k, \lambda^k)$ it sets

$$\mathcal{A}_k = \{ i : \lambda^k_i + c(y^k_i - \psi_i) > 0 \}, \qquad \mathcal{I}_k = \{ i : \lambda^k_i + c(y^k_i - \psi_i) \le 0 \},$$

and determines $(y^{k+1}, \lambda^{k+1})$ from

$$Ay^{k+1} + \lambda^{k+1} = f, \qquad y^{k+1} = \psi \text{ on } \mathcal{A}_k, \qquad \lambda^{k+1} = 0 \text{ on } \mathcal{I}_k .$$

For a rectangular matrix $B$, $B_+$ is the matrix of the positive parts $\max(0, b_{ij})$ of its entries, and $\|\cdot\|_1$ is the matrix norm subordinate to the one-norms, i.e. the maximal absolute column sum. Thus $\|B_+\|_1 = \max_j \sum_i \max(0, b_{ij})$.

Formalization targets

Goal: Theorem 3.3

If $A$ is a P-matrix and for every partition of $\{1, \dots, n\}$ into disjoint sets $\mathcal{I}$ and $\mathcal{A}$

$$\big\|(A_{\mathcal{I}}^{-1} A_{\mathcal{I}\mathcal{A}})_+\big\|_1 < 1 \qquad\text{and}\qquad \sum_{i \in \mathcal{I}} (A_{\mathcal{I}}^{-1} v)_i \ge 0 \ \text{ for all } v \ge 0,$$

then for every solution $x^* = (y^*, \lambda^*)$ of (3.1) and every run of the algorithm, $x^k = (y^k, \lambda^k) \to x^*$.

Milestones

The proof in the paper passes through the following claims, in order: the sign facts $\lambda^k_{\mathcal{I}_k} \le 0$ and $y^k_{\mathcal{A}_k} \ge \psi_{\mathcal{A}_k}$ for $k \ge 1$; the update formula on the inactive set derived from the Newton system (3.3); inequality (3.6); inequality (3.7),

$$\sum_{i=1}^n (y^{k+1}_i - y^k_i) \le -\|y^k - \psi\|_{1,\mathcal{A}_k} + \big\|(A_{\mathcal{I}_k}^{-1} A_{\mathcal{I}_k\mathcal{A}_k})_+\big\|_1 \, \|y^k - \psi\|_{1,\mathcal{A}_k};$$

and the fact that a repeated active set yields a solution of (3.1).

Companion statements

Existence of a run from every initial pair when $A$ is a P-matrix; Remark 3.3 (an M-matrix satisfies the conditions of Theorem 3.3); and the stronger conclusion of the proof, that every run reaches $x^*$ after finitely many steps.

Significance

Theorem 3.3 gives global convergence of the primal-dual active set strategy for a class of matrices that is not restricted to M-matrices, and identifies $\mathcal{M}(y) = \sum_i y_i$ as a merit function under its two conditions. The paper points out that the result applies to discretizations of control-constrained optimal control problems, where the system matrix is generally not an M-matrix, and Remark 3.3 shows that the M-matrix case falls under it.

The result is proved on paper; none of it has been machine-checked. A formal development would produce a checked global convergence theorem for a semismooth Newton method, together with reusable statements about principal blocks of P-matrices and the subordinate one-norm of a positive part. The step from "the merit function does not increase" to convergence is the part of the paper's argument that a formal proof has to make fully precise (see below).

Difficulty

The obvious route is the paper's: show that $\sum_i y^k_i$ decreases strictly until the active set repeats, conclude that no active set can recur, and use finiteness of the set of possible active sets. The inequalities (3.6) and (3.7) give only a non-strict decrease. The paper's assertion "$< 0$, unless $y^{k+1} = y^k$" does not follow from them as written: the right-hand side of (3.7) vanishes when $y^k = \psi$ on $\mathcal{A}_k$, and the cone condition permits $\sum_{i \in \mathcal{I}} (A_{\mathcal{I}}^{-1}\lambda^k_{\mathcal{I}})_i = 0$ with $\lambda^k_{\mathcal{I}} \neq 0$. A complete proof has to exclude cycling among active sets at equal merit value by a further argument. That is the central difficulty of this mission, and the reason the strict part of (3.7) is not stated as a milestone.

Formalization scope

Vectors are Fin n → ℝ and matrices Matrix (Fin n) (Fin n) ℝ; index sets are Finset (Fin n) and blocks are submatrices along the subtype inclusion. The P-matrix predicate is the published definition RobinsonSR.Schur.IsPMatrix (every nonempty principal minor is positive). The algorithm is encoded as a relation between consecutive iterates, and a run is an infinite sequence related step by step with arbitrary initial data; the stopping option of the algorithm is not modelled. The solution of (3.1) enters as a hypothesis: for a P-matrix it exists and is unique, so no assumption beyond the paper is made. Convergence is convergence of $(y^k, \lambda^k)$ in $\mathbb{R}^n \times \mathbb{R}^n$. The parameter $c > 0$ is the paper's standing assumption.

Conventions committed to:

• both conditions of Theorem 3.3 range over every index set, including $\mathcal{I} = \emptyset$ and $\mathcal{I} = \{1,\dots,n\}$;
• $\|B_+\|_1$ is the maximal column sum of the positive part, computed as a supremum in $\mathbb{R}_{\ge 0}$, and equals $0$ for a matrix with no columns;
• the cone condition is $\sum_{i\in\mathcal{I}}(A_{\mathcal{I}}^{-1}v)_i \ge 0$ for $v \ge 0$, not entrywise nonnegativity of $A_{\mathcal{I}}^{-1}$;
• the milestones on the iterates hold from $k = 1$ on, as in the paper.

The goal does not assume a decrease of the merit function or finite termination; a formalization that adds either as a hypothesis would be trivial and is ruled out. The run-existence companion shows that the hypothesis "let $(y^k,\lambda^k)$ be a run" can be met, so the goal is not vacuous.

A complete development needs: invertibility of principal blocks of a P-matrix, the block form of one step, the elementary estimate of $\sum_i (B x)_i$ by $\|B_+\|_1 \sum_j x_j$ for $x \ge 0$, and an argument that excludes cycling. Proofs of any milestone, of the companion statements, or of the goal through an argument other than the paper's are welcome.

Selected references

• M. Hintermüller, K. Ito, K. Kunisch, The primal-dual active set strategy as a semismooth Newton method, SIAM J. Optim. 13(3):865–888, 2002. https://doi.org/10.1137/S1052623401383558 ; authors' version https://hal.science/hal-01660511v1
• A. Berman, R. J. Plemmons, Nonnegative Matrices in the Mathematical Sciences, Academic Press, 1979 (SIAM Classics reprint 1994). https://doi.org/10.1137/1.9781611971262

Motivation

Optimal control problems with pointwise bounds on the control, such as

$$\min_u \tfrac12\|y - z\|_{L^2}^2 + \tfrac{\beta}{2}\|u\|_{L^2}^2 \quad\text{subject to } -\Delta y = u \text{ in } \Omega,\ y = 0 \text{ on } \partial\Omega,\ u \le \psi,$$

reduce, after eliminating the state, to a quadratic program in $L^2(\Omega)$ with a pointwise upper bound. The primal-dual active set strategy solves such problems by guessing the set where the bound is active, solving one linear system on that guess, and updating the guess from the new primal and dual values. It is used routinely in PDE-constrained optimization because each iteration costs one linear solve and, in practice, very few iterations are needed.

Hintermüller, Ito and Kunisch (SIAM J. Optim. 13 (2002) 865–888; authors' version hal-01660511) showed that this strategy is a semismooth Newton method for a nonsmooth reformulation of the optimality system, built on the notion of slant differentiability of Chen, Nashed and Qi (SIAM J. Numer. Anal. 38 (2000) 1200–1216). In finite dimension this gives local superlinear convergence (missions I–IV of this series). This mission concerns the paper's function-space result, Theorem 4.1: in $L^2(\Omega)$, under a norm gap on the operator, the method still converges locally superlinearly. Ulbrich (SIAM J. Optim. 13 (2003) 805–841) developed a related semismooth Newton theory for superposition operators in function spaces, which also requires a norm gap.

Setting

Let $(\Omega, \mu)$ be a finite measure space (the paper takes a bounded domain $\Omega \subset \mathbb{R}^n$ with Lebesgue measure) and $L^p = L^p(\Omega)$ the real Lebesgue spaces, with inner product $(\cdot,\cdot)$ on $L^2$. The paper writes $A$ for an operator and calligraphic $\mathcal{A}$, $\mathcal{I}$ for active and inactive sets.

Slant differentiability (Definition 1). For Banach spaces $X$, $Z$, a map $F : X \to Z$ is slantly differentiable in an open set $U$ if there is $G : U \to \mathcal{L}(X, Z)$, a slanting function, with

$$\lim_{h\to 0}\frac{1}{\|h\|}\,\|F(x+h) - F(x) - G(x+h)h\| = 0 \qquad (x \in U). \tag{A}$$

For $\delta \in \mathbb{R}$, the candidate slanting function of the pointwise max is $G_m(y)(x) = 1$ if $y(x) > 0$, $0$ if $y(x) < 0$, $\delta$ if $y(x) = 0$ (4.1).

The problem. Given $f, \psi \in L^2$ and a self-adjoint $A \in \mathcal{L}(L^2)$ with

$$(Ay, y) \ge \gamma\|y\|^2 \quad (\gamma > 0), \tag{H1}$$

minimize $\tfrac12(y, Ay) - (f, y)$ subject to $y \le \psi$. The solution $y^*$ and its multiplier $\lambda^*$ are the unique solution of

$$Ay^* + \lambda^* = f, \qquad \lambda^* - \max(0, \lambda^* + c(y^* - \psi)) = 0 \ \text{a.e.}, \tag{4.2}$$

for any fixed $c > 0$. The norm gap hypothesis is

$$A = C + \beta I, \qquad C \in \mathcal{L}(L^2, L^q),\quad \beta > 0,\quad q > 2. \tag{H2}$$

The algorithm. From $(y^0, \lambda^0) \in L^2 \times L^2$, set $\mathcal{A}_k = \{x : \lambda^k(x) + c(y^k(x) - \psi(x)) > 0\}$, $\mathcal{I}_k = \Omega\setminus\mathcal{A}_k$, and solve

$$Ay^{k+1} + \lambda^{k+1} = f, \qquad y^{k+1} = \psi \text{ on } \mathcal{A}_k, \qquad \lambda^{k+1} = 0 \text{ on } \mathcal{I}_k.$$

A sequence obeying this rule for every $k$ is a run. Superlinear convergence of $x^k \to x^*$ means $x^k \to x^*$ and, for every $\eta > 0$, eventually $\|x^{k+1} - x^*\| \le \eta\|x^k - x^*\|$.

Formalization targets

Goal: Theorem 4.1 (p. 15)

The page reads: "Assume that (H1), (H2) hold and that ψ and f are in L^q(Ω). Then the primal-dual active set strategy or, equivalently, the semismooth Newton method converge superlinearly if ‖y⁰ − y*‖ is sufficiently small and λ⁰ = β(y⁰ − ψ)." The formalized goal is the statement the proof establishes: with $c = \beta$, there is $\rho > 0$ such that every run with

$$\|y^0 - y^*\|_{L^2} < \rho, \qquad \lambda^0 = f - Ay^0$$

converges superlinearly to $(y^*, \lambda^*)$ in $L^2 \times L^2$. The initialization differs from the printed one; see Formalization scope.

Milestones

1. Theorem 1.1 (Chen–Nashed–Qi): slant differentiability near a zero $x^*$ with uniformly bounded $G(x)^{-1}$ gives local superlinear convergence of $x^{k+1} = x^k - G(x^k)^{-1}F(x^k)$.
2. Proposition 4.1(ii): for $1 \le p < q \le \infty$, $\max(0,\cdot) : L^q \to L^p$ is slantly differentiable with slanting function $G_m$.
3. (4.9)–(4.10): for $c = \beta$, (4.2) is equivalent to $\beta y^* - \beta\psi + \max(0, Cy^* - f + \beta\psi) = 0$ and $\lambda^* = f - Cy^* - \beta y^*$.
4. Equivalence of the algorithms: along a run with $c = \beta$ and $\lambda^0 = f - Ay^0$, $\lambda^k + \beta(y^k - \psi) = f - Cy^k - \beta\psi$ for all $k$, so the active sets are those of the reduced algorithm.
5. Slant differentiability of the reduced map $F(y) = \beta y - \beta\psi + \max(0, Cy - f + \beta\psi)$ on $L^2$, with $G_F(y+h) = \beta I + G_m(Cy - f + \beta\psi + Ch)C$, $\delta = 1$.
6. Bounded inverses: $G_F(z)$ is invertible on $L^2$, with $\|G_F(z)^{-1}\|$ bounded uniformly in $z$.

Two companion items accompany the goal: existence of a run from every $y^0$, and Proposition 4.1(i), the paper's counterexamples showing that $G_m$ fails to be a slanting function from $L^p$ to $L^p$.

Significance

Theorem 4.1 is the infinite-dimensional justification of a method used throughout control-constrained optimal control: it shows that the method converges locally superlinearly for the continuous problem itself, not only for each discretization. Example 1 of the paper shows that distributed control of the Poisson equation, and boundary control with Sobolev embeddings, satisfy (H2). Proposition 4.1 isolates the structural reason: the pointwise max is a Newton-differentiable map only between spaces with a norm gap.

None of these results is known to be machine-checked. A complete development would provide slant differentiability in Banach spaces, the Chen–Nashed–Qi local convergence theorem, the norm-gap differentiability of $\max(0,\cdot)$ between Lebesgue spaces, and invertibility of block operators on $L^2$ of a measurable partition. Each of these is reusable for other semismooth Newton analyses in function space.

Difficulty

The obvious argument fails at the max. In finite dimension, $\max(0,\cdot)$ is slantly differentiable with $G_m$ because a small perturbation does not change the sign of any nonzero component. In $L^2$, a perturbation $h$ of small norm can flip the sign of $y$ on a set of positive measure, and on that set the error $\max(0, y+h) - \max(0, y) - G_m(y+h)h$ is as large as $|y|$. Proposition 4.1(i) shows that this error is not $o(\|h\|_{L^p})$ in $L^p$ for any $p$. Superlinear convergence in $L^2$ therefore needs the argument of the max to be controlled in a stronger norm than the error is measured in, which is what (H2) and $f, \psi \in L^q$ provide. A second difficulty is that the two-variable algorithm is not literally a Newton method for a map of $y$ alone; it must be identified with the reduced algorithm, and that identification depends on the initialization.

Formalization scope

• Spaces. $L^p$ is Mathlib's Lp ℝ p μ on an arbitrary finite measure space; $q \in (2, \infty]$ is an extended real with 2 < q. The paper's bounded Lipschitz domain is replaced by a finite measure space: the proofs use only $|\Omega| < \infty$, and Mathlib has no notion of Lipschitz boundary.
• Pointwise statements on $L^p$ classes are almost everywhere statements about representatives. (H2) is stated as $Ay = Cy + \beta y$ a.e. for every $y \in L^2$, which avoids building the inclusion $L^q \subseteq L^2$.
• Hypotheses kept: self-adjointness of $A$, (H1) with explicit $\gamma > 0$, (H2) with $\beta > 0$ and $q > 2$, $f, \psi \in L^q$. The solution $(y^*, \lambda^*)$ of (4.2) is a hypothesis; it exists and is unique under (H1), as the paper states on p. 12.
• The initialization. The printed $\lambda^0 = \beta(y^0 - \psi)$ does not support the proof, which uses "$\lambda^0 + \beta(y^0 - \psi) = f - Cy^0 - \beta\psi$" and $c = \beta$. That identity holds exactly for $\lambda^0 = f - Ay^0$. The goal uses $\lambda^0 = f - Ay^0$ and $c = \beta$. With the printed initialization the first active set is $\{y^0 > \psi\}$, and for $y^0$ slightly below $y^*$ the first iterate solves the unconstrained problem, which can be far from $y^*$.
• Superlinear convergence is stated for the pair $(y^k, \lambda^k)$ in the product norm $\max(\|y\|, \|\lambda\|)$, without dividing by $\|x^k - x^*\|$, so finite termination qualifies. "Sufficiently small" is an $\exists\rho > 0$ chosen before the run, on $\|y^0 - y^*\|$ only.
• No trivialization. The goal mentions neither the reduced algorithm, nor $G_F$, nor Proposition 4.1; the norm gap $q > 2$ is not weakened to $C \in \mathcal{L}(L^2)$, under which Proposition 4.1(i) shows that $G_m$ need not be a slanting function; and the companion run_exists shows that the runs the goal quantifies over exist.
• Definitions (PDASNewton.FunSpace.Setting): slanting functions, superlinear convergence, $G_m$, property (A) for $\max(0,\cdot) : L^q \to L^p$, system (4.2), one step and runs of the algorithm. These duplicate the finite-dimensional definitions of the other missions of the series in a separate namespace.

Proofs of any milestone are welcome, as are reusable lemmas on Lebesgue spaces (Hölder on sets of small measure, multiplication operators, Lax–Milgram on $L^2$ of a measurable subset).

Selected references

• M. Hintermüller, K. Ito, K. Kunisch, The primal-dual active set strategy as a semismooth Newton method, SIAM J. Optim. 13(3) (2002) 865–888. https://doi.org/10.1137/S1052623401383558 ; authors' version https://hal.science/hal-01660511
• X. Chen, Z. Nashed, L. Qi, Smoothing methods and semismooth methods for nondifferentiable operator equations, SIAM J. Numer. Anal. 38 (2000) 1200–1216. https://doi.org/10.1137/S0036142999356719
• M. Ulbrich, Semismooth Newton methods for operator equations in function spaces, SIAM J. Optim. 13 (2003) 805–841. https://doi.org/10.1137/S1052623400371569
• L. Qi, J. Sun, A nonsmooth version of Newton's method, Math. Programming 58 (1993) 353–367. https://doi.org/10.1007/BF01581275

Motivation

Many problems in optimization and in the discretization of partial differential equations reduce to a linear complementarity problem with an obstacle: find $y, \lambda \in \mathbb{R}^n$ with

$$Ay + \lambda = f, \qquad y \le \psi, \quad \lambda \ge 0, \quad (\lambda, y - \psi) = 0,$$

for a matrix $A$ and vectors $f, \psi$. When $A$ is symmetric positive definite this is the optimality system of the quadratic program $\min \tfrac12 (y, Ay) - (f, y)$ subject to $y \le \psi$. Finite-difference and finite-element discretizations of obstacle problems and of control-constrained optimal control problems lead to systems of exactly this form, with $n$ in the thousands or millions.

The primal-dual active set strategy solves such systems by guessing, at every iteration, which constraints are active and solving one linear system for that guess. It was introduced as a special case of generalized Moreau–Yosida approximations by Ito and Kunisch, and was observed to converge in a few iterations, typically without a line search. M. Hintermüller, K. Ito and K. Kunisch, The primal-dual active set strategy as a semismooth Newton method, SIAM J. Optim. 13(3):865–888, 2002 (DOI 10.1137/S1052623401383558; authors' version HAL hal-01660511) explained that behaviour: the strategy is a semismooth Newton method for a nonsmooth reformulation of the problem, and therefore inherits local superlinear convergence. This mission formalizes that local theory, the first of a five-part series on the paper.

Setting

Fix $n \ge 0$, $A \in \mathbb{R}^{n\times n}$, $f, \psi \in \mathbb{R}^n$ and $c > 0$. Vectors are ordered componentwise. The complementarity conditions are equivalent to one nonsmooth equation, so the problem becomes

$$Ay + \lambda = f, \qquad \lambda - \max(0, \lambda + c(y - \psi)) = 0, \tag{3.1}$$

with the maximum taken componentwise. The paper writes $A$ for the matrix and calligraphic $\mathcal{A}$, $\mathcal{I}$ for index sets; the two are unrelated.

Primal-dual active set algorithm. Given an iterate $(y^k, \lambda^k)$, the active set is $\mathcal{A}_k = \{i : \lambda^k_i + c(y^k - \psi)_i > 0\}$ and the inactive set is $\mathcal{I}_k = \{i : \lambda^k_i + c(y^k - \psi)_i \le 0\}$. The next iterate solves

$$Ay^{k+1} + \lambda^{k+1} = f, \qquad y^{k+1} = \psi \text{ on } \mathcal{A}_k, \qquad \lambda^{k+1} = 0 \text{ on } \mathcal{I}_k.$$

A run is a sequence of iterates, from arbitrary initial data $(y^0, \lambda^0)$, in which every pair follows its predecessor in this way.

A P-matrix is a square matrix all of whose principal minors are positive. For a P-matrix every principal block $A_{\mathcal{I}}$ is invertible, so every step is uniquely solvable, and (3.1) has exactly one solution $(y^*, \lambda^*)$.

Slant differentiability. Let $X, Z$ be Banach spaces and $F : X \to Z$. A map $G$ from $X$ to the bounded linear operators $\mathcal{L}(X, Z)$ is a slanting function for $F$ in an open set $U$ if for every $x \in U$

$$\lim_{h \to 0} \frac{1}{\|h\|}\,\|F(x+h) - F(x) - G(x+h)h\| = 0. \tag{A}$$

The derivative is evaluated at the perturbed point $x + h$, not at $x$; this is what lets a nonsmooth map such as $\max(0,\cdot)$ qualify. A sequence $x^k$ converges superlinearly to $x^*$ if $x^k \to x^*$ and for every $\eta > 0$, eventually $\|x^{k+1} - x^*\| \le \eta\,\|x^k - x^*\|$.

For (3.1) put $F(y, \lambda) = \bigl(Ay + \lambda - f,\ \lambda - \max(0, \lambda + c(y-\psi))\bigr)$. For $\delta \in \mathbb{R}^n$ let $G_m(y)$ be the diagonal matrix with entries $0$ where $y_i < 0$, $1$ where $y_i > 0$ and $\delta_i$ where $y_i = 0$. The paper's $G_F(x)$ is the block matrix (2.4) obtained by using $G_m$ with $\delta = 0$ for the max-term of $F$.

Formalization targets

Goal: Theorem 3.1

For a P-matrix $A$, $c > 0$ and the solution $x^* = (y^*, \lambda^*)$ of (3.1), there is $\rho > 0$ such that every run with $\|x^0 - x^*\| < \rho$ satisfies

$$x^k \to x^*, \qquad \forall \eta > 0:\ \|x^{k+1} - x^*\| \le \eta\,\|x^k - x^*\| \text{ for all large } k.$$

The radius depends on the data and the solution, not on the run. The goal speaks about the active set algorithm itself; its identification with a semismooth Newton method is a milestone.

Milestones

1. Theorem 1.1 (cited from Chen, Nashed and Qi): if $F(x^*) = 0$, $G$ is a slanting function for $F$ on an open neighbourhood $U$ of $x^*$, and $G(x)$ is invertible for $x \in U$ with $\{\|G(x)^{-1}\|\}$ bounded, then the Newton iteration $x^{k+1} = x^k - G(x^k)^{-1}F(x^k)$ converges superlinearly to $x^*$ from every $x^0$ close enough to $x^*$.
2. Lemma 3.1: $y \mapsto \max(0, y)$ on $\mathbb{R}^n$ is slantly differentiable with slanting function $G_m$, for every $\delta$.
3. §3, p. 7: $F$ is slantly differentiable on $\mathbb{R}^n \times \mathbb{R}^n$ with slanting function $G_F$.
4. §2, (2.7)–(2.8): for a P-matrix, every $G_F(x)$ is invertible and the inverses are uniformly bounded.
5. §2, (2.4)–(2.6): $G_F(x^k)(x^{k+1} - x^k) = -F(x^k)$ holds exactly when $x^{k+1}$ follows $x^k$ by one active set step.

Three companion statements accompany the goal: a run exists from every initial pair; a run that converges to $x^*$ reaches it after finitely many steps (§3, p. 7); and if $y^* < \psi$ and $\lambda^0 + c(y^0 - \psi) \le 0$, the run equals $x^*$ from the first iterate on (Remark 3.4 (iii)).

Significance

The result. Theorem 3.1 explains why the primal-dual active set strategy converges fast in practice, and it does so through a notion, slant differentiability, that carries over to function spaces, where Clarke's generalized Jacobian is not available. The paper's later sections (missions II–V of this series) use the same framework for global convergence under M-matrix, diagonal-dominance and perturbation hypotheses, and for superlinear convergence in $L^2(\Omega)$. Theorem 1.1 itself is an abstract local convergence result for nonsmooth Newton methods that applies far beyond complementarity problems.

Formalizing it. The results are proved on paper; none of them is machine-checked. Mathlib has the Fréchet derivative and Newton-type fixed-point arguments, but no slant differentiability, no semismooth Newton method and no theory of the primal-dual active set strategy. A completed mission gives a reusable Lean statement and proof of the Chen–Nashed–Qi theorem in Banach spaces, slant differentiability of the componentwise maximum, and a verified bridge between a combinatorial active set iteration and a Newton iteration.

Difficulty

The map $\max(0, \cdot)$ is not differentiable at $0$, so the classical Newton convergence theorem does not apply to $F$. The obvious substitute is a condition with the derivative at the base point, $\|F(x) - F(x-h) - G(x)h\| = o(\|h\|)$; the paper notes (p. 3) that with it Theorem 1.1 would need an additional uniformity assumption in $x \in U$ when $X$ is infinite dimensional. Statements and proofs therefore have to work with the derivative at the perturbed point, which is the point at which the Newton iteration evaluates it. A second difficulty is bookkeeping: the identification of the Newton step with the active set step is a block computation over index sets that change with $k$, and the bound on $\|G_F(x)^{-1}\|$ must hold uniformly over all $2^n$ active sets, using only that principal blocks of a P-matrix are invertible.

Formalization scope

Vectors are Fin n → ℝ, matrices Matrix (Fin n) (Fin n) ℝ, and the order is componentwise. Pairs $(y, \lambda)$ live in (Fin n → ℝ) × (Fin n → ℝ) with the max of the two sup norms; in finite dimension every norm gives the same notion of superlinear convergence. Property (A) is stated as a little-o relation, so nothing is divided by $\|h\|$; boundedness of $\{G(x)\}$, required in the original Chen–Nashed–Qi definition, is deliberately not part of it, as in the paper. $F$ and $G$ are total maps; the paper's domain $D \supseteq U$ plays no role.

Standing assumptions: $c > 0$ in every statement; $A$ is a P-matrix wherever the paper uses it (the goal, the inverse bound, the existence of runs, finite termination, Remark 3.4 (iii)), via the published definition RobinsonSR.Schur.IsPMatrix. The solution $(y^*, \lambda^*)$ of (3.1) is a hypothesis; for a P-matrix it exists and is unique, which is a separate classical theorem not restated here. A step of the algorithm is a relation and a run is an infinite sequence of steps; the "Stop" option of the algorithm is not modelled. Nonsingularity in Theorem 1.1 means a two-sided bounded linear inverse, and the Newton iteration is written as the linear equation $G(x^k)(x^{k+1} - x^k) = -F(x^k)$.

The radius $\rho$ in Theorem 1.1 and Theorem 3.1 is quantified before the sequence; a formalization in which $\rho$ may depend on the run, or in which superlinear convergence drops either the convergence or the rate condition, does not state the theorem. A companion item shows runs exist, so the goal is not vacuous.

Contributions welcome: proofs of Theorem 1.1 in Banach spaces (reusable for missions II–V), of Lemma 3.1, of the block identities, and of the uniform inverse bound over all active sets.

Selected references

• M. Hintermüller, K. Ito, K. Kunisch, The primal-dual active set strategy as a semismooth Newton method, SIAM J. Optim. 13(3):865–888, 2002. doi:10.1137/S1052623401383558; authors' version hal-01660511v1.
• X. Chen, Z. Nashed, L. Qi, Smoothing methods and semismooth methods for nondifferentiable operator equations, SIAM J. Numer. Anal. 38(4):1200–1216, 2000. doi:10.1137/S0036142999356719
• L. Qi, J. Sun, A nonsmooth version of Newton's method, Math. Program. 58:353–367, 1993. doi:10.1007/BF01581275
• A. Berman, R. J. Plemmons, Nonnegative Matrices in the Mathematical Sciences, SIAM Classics in Applied Mathematics 9, 1994. doi:10.1137/1.9781611971262

Motivation

Many problems in signal processing, statistics and decision theory have two blocks of variables coupled by a smooth term: matrix factorization, blind deconvolution, feasibility problems for two sets, and models of two agents adjusting their decisions in turn. Minimizing over one block at a time (alternating minimization, or Gauss–Seidel) is the natural method. Without convexity it can oscillate, and its iterates need not converge. Attouch, Bolte, Redont and Soubeyran (Math. Oper. Res. 35(2), 2010; arXiv:0801.1780) regularize each block step by a proximal term. They prove that, when the objective satisfies the Kurdyka–Łojasiewicz (KL) inequality, every sequence of the method that does not run off to infinity has finite length and converges to a critical point.

Timeline:

• 1963: Łojasiewicz proves the gradient inequality for real-analytic functions, which gives convergence of bounded gradient trajectories.
• 1998: Kurdyka extends the inequality to differentiable functions definable in an o-minimal structure.
• 2006–2007: Bolte, Daniilidis, Lewis and Shiota extend it to nonsmooth lower semicontinuous subanalytic and definable functions, with the limiting subdifferential.
• 2009: Attouch and Bolte use the Łojasiewicz inequality to prove convergence of the proximal point method for nonconvex functions (Math. Program. 116, 2009).
• 2010: this paper treats the alternating proximal scheme with a general smooth coupling $Q$.
• 2014: Bolte, Sabach and Teboulle replace the exact block minimizations by linearized proximal steps (PALM, Math. Program. 146, 2014).

Setting

Let $f:\mathbb R^n\to\mathbb R\cup\{+\infty\}$ and $g:\mathbb R^m\to\mathbb R\cup\{+\infty\}$ be proper (never $-\infty$, finite somewhere) and lower semicontinuous. Let $Q:\mathbb R^n\times\mathbb R^m\to\mathbb R$ be continuously differentiable, with $\nabla Q$ Lipschitz on bounded sets. This is assumption (H) for

$$L(x,y)=f(x)+Q(x,y)+g(y).$$

The product $\mathbb R^n\times\mathbb R^m$ carries the Euclidean norm $\|(x,y)\|^2=\|x\|^2+\|y\|^2$. Write $\nabla Q=(\nabla_xQ,\nabla_yQ)$.

Given $(x_0,y_0)$ and step sizes $\lambda_k,\mu_k>0$, the proximal alternating scheme is

$$x_{k+1}\in\operatorname{argmin}_u\Big\{L(u,y_k)+\tfrac1{2\lambda_k}\|u-x_k\|^2\Big\},\qquad y_{k+1}\in\operatorname{argmin}_v\Big\{L(x_{k+1},v)+\tfrac1{2\mu_k}\|v-y_k\|^2\Big\}.\tag{5–6}$$

Assumption (H1) asks that $\inf L>-\infty$, that $L(\cdot,y_0)$ be proper, and that $\lambda_k,\mu_k\in(r_-,r_+)$ for all $k$, for some $0<r_-<r_+$.

The Fréchet subdifferential $\hat\partial L(z)$ is the set of $v$ with $\liminf_{w\to z,\,w\ne z}\big(L(w)-L(z)-\langle v,w-z\rangle\big)/\|w-z\|\ge0$. The limiting subdifferential $\partial L(z)$ collects the limits of $v_j\in\hat\partial L(z_j)$ along $z_j\to z$ with $L(z_j)\to L(z)$. A point $z$ is critical if $0\in\partial L(z)$. $L$ has the KL property at $\bar z\in\operatorname{dom}\partial L$ if there are $\eta>0$, a neighbourhood $U$ of $\bar z$ and a continuous concave $\varphi:[0,\eta)\to\mathbb R_+$ with $\varphi(0)=0$, $\varphi\in C^1(0,\eta)$, $\varphi'>0$, such that

$$\varphi'\big(L(z)-L(\bar z)\big)\,\operatorname{dist}\big(0,\partial L(z)\big)\ge1\quad\text{whenever }z\in U,\ L(\bar z)<L(z)<L(\bar z)+\eta .$$

Formalization targets

Goal: Theorem 9 (p. 11)

If $L$ satisfies (H), (H1) and has the KL property at every point of $\operatorname{dom}\partial L$, then every sequence of (5–6) satisfies

$$\|(x_k,y_k)\|\to\infty\qquad\text{or}\qquad\Big(\sum_k\|x_{k+1}-x_k\|+\|y_{k+1}-y_k\|<\infty\ \text{ and }\ (x_k,y_k)\to z^\star,\ 0\in\partial L(z^\star)\Big).$$

Milestones

• Proposition 3 (p. 5): $\partial L(x,y)=\{\partial f(x)+\nabla_xQ(x,y)\}\times\{\partial g(y)+\nabla_yQ(x,y)\}=\partial_xL(x,y)\times\partial_yL(x,y)$ on $\operatorname{dom}L$.
• Lemma 5 (p. 6), in four parts: the scheme is well defined; the decrease estimate (7); summability of the squared steps; the explicit subgradient $(x_k^*,y_k^*)\in\partial L(x_k,y_k)$ (8), which tends to $0$ along bounded subsequences.
• Proposition 6 (p. 7): limit points are critical, and $L$ is finite and constant on them.
• (19)–(20), (24), (25), (26) from the proof of Theorem 8 (pp. 9–11).
• Theorem 8 (p. 9): a local convergence result with explicit constants: $M=2r_+(C+1/r_-)$, the trapping radius (16), and the tail bound (18).

Significance

Theorem 9 makes the KL inequality the single hypothesis behind global convergence of a nonconvex block method. Semi-algebraic, subanalytic and definable data satisfy it, so the theorem covers sparse and low-rank regularizers, indicator functions of semi-algebraic sets, and polynomial couplings. With $f=\delta_C$, $g=\delta_D$ and $Q(x,y)=\tfrac12\|x-y\|^2$ it yields convergence of an averaged alternating projection method between two closed sets (§3.4 of the paper). The proof pattern of sufficient decrease, a relative-error subgradient bound and the KL inequality was later abstracted for PALM and many other algorithms. A formal proof of this paper's argument therefore covers the template those analyses reuse.

The result has been proved since 2010. To our knowledge it has not been formalized in any proof assistant. Mathlib has the Fréchet derivative and the gradient, but no limiting subdifferential and no KL theory. The published Prove2Me definitions of the limiting subdifferential and of the KL property (from the Li–Pong ADMM series) are reused here, so this mission also tests that infrastructure on a second algorithm.

Difficulty

The decrease estimate (7) gives only $\sum\|z_{k+1}-z_k\|^2<\infty$. That makes steps vanish and limit points critical (Proposition 6), but it does not give convergence, because square-summable steps can have infinite total length. Turning square summability into summability needs the KL inequality. That inequality holds only near a critical point and only on a sublevel band, so the iterates must first be shown never to leave a ball $B(\bar z,\rho)$, by an induction that uses the very finite-length estimate it is trying to establish. The constants in (16) are arranged so that this induction closes. A second difficulty is nonsmooth calculus: the subdifferential sum rule for $f(x)+g(y)+Q(x,y)$ and the closedness of the graph of $\partial L$ along sequences with converging function values. Lower semicontinuity alone does not give convergence of the function values; the proximal inequalities do.

Formalization scope

Everything lives in the namespace ProxAltMin.Conv, in Lean 4 with Mathlib. Conventions:

• $\mathbb R^n$, $\mathbb R^m$ are EuclideanSpace ℝ (Fin n) and EuclideanSpace ℝ (Fin m), and the product is the $L^2$ product WithLp 2 (… × …), not the sup-norm product.
• $\mathbb R\cup\{+\infty\}$ is EReal with "never $-\infty$" as part of properness.
• The argmin steps are inequalities quantified over all competitors, so every selection of minimizers is covered.
• Balls are open.
• The KL data $(\eta,U,\varphi)$ of Theorem 8 are explicit parameters, with real $\eta>0$; this is equivalent to $\eta\in(0,+\infty]$.
• Inside the proof of Theorem 8 the paper normalizes $L(\bar z)=0$; the milestones taken from that proof carry $l_i-\bar l$ instead.
• Every infinite sum in a conclusion comes with an explicit summability claim.
• Theorem 9's "KL property at each point of the domain of $f$" is read as "at each point of $\operatorname{dom}\partial L$", the domain on which Definition 7 places the KL property.

Theorem and page numbers are those of arXiv:0801.1780v3 (22 Jan 2013), which is the version read; the published article has the same content.

A trivializing formalization is ruled out: the run predicate is satisfiable (Lemma 5's existence milestone; for $f=g=Q=0$ any constant sequence is a run), the constants $M$ and $C$ are tied to the data rather than chosen by the prover, and no conclusion is a sum or a distance whose library value is $0$ on bad input.

Needed infrastructure: the subdifferential calculus behind Proposition 3 (a smooth perturbation and a separable sum), the closedness of the limiting subdifferential, existence of proximal minimizers for coercive lower semicontinuous functions on $\mathbb R^n$, and the elementary but delicate bookkeeping of Theorem 8. The first two are reusable well beyond this mission. Contributions welcome: proofs of any milestone, and general lemmas on limiting subdifferentials stated as separate theorems.

Selected references

• H. Attouch, J. Bolte, P. Redont, A. Soubeyran, Proximal alternating minimization and projection methods for nonconvex problems: an approach based on the Kurdyka–Łojasiewicz inequality, Mathematics of Operations Research 35(2):438–457, 2010. https://doi.org/10.1287/moor.1100.0449 (arXiv:0801.1780v3, https://arxiv.org/abs/0801.1780v3)
• H. Attouch, J. Bolte, On the convergence of the proximal algorithm for nonsmooth functions involving analytic features, Mathematical Programming 116:5–16, 2009. https://doi.org/10.1007/s10107-007-0133-5
• J. Bolte, A. Daniilidis, A. Lewis, The Łojasiewicz inequality for nonsmooth subanalytic functions with applications to subgradient dynamical systems, SIAM J. Optim. 17(4):1205–1223, 2007. https://doi.org/10.1137/050644641
• K. Kurdyka, On gradients of functions definable in o-minimal structures, Annales de l'Institut Fourier 48(3):769–783, 1998. https://doi.org/10.5802/aif.1638
• J. Bolte, S. Sabach, M. Teboulle, Proximal alternating linearized minimization for nonconvex and nonsmooth problems, Mathematical Programming 146:459–494, 2014. https://doi.org/10.1007/s10107-013-0701-9
• R. T. Rockafellar, R. J.-B. Wets, Variational Analysis, Springer, 1998. https://doi.org/10.1007/978-3-642-02431-3

Motivation

Many problems in statistics, signal processing and machine learning minimize the sum of a smooth function and a convex but nonsmooth one: $\ell_1$-regularized least squares (the Lasso), bound-constrained problems, and group-sparse regression are standard examples. Tseng and Yun (Math. Program. Ser. B 117 (2009) 387–423) proposed the coordinate gradient descent (CGD) method for this class and proved that it converges linearly under a local Lipschitzian error bound, Assumption 2(a) of their paper. An assumption is only as useful as the list of problems known to satisfy it. Section 6 of the paper supplies that list, and this mission formalizes it.

Error bounds of this kind have a history in smooth constrained optimization. Luo and Tseng proved them for quadratic objectives over polyhedral sets (Luo–Tseng 1992, SIAM J. Optim.), for strongly convex functions composed with a linear map (Luo–Tseng 1992, SIAM J. Control Optim.), and for dual functionals of linearly constrained strictly convex programs (Luo–Tseng 1993, Math. Oper. Res.), and used them to derive linear rates for feasible descent methods (Luo–Tseng 1993, Ann. Oper. Res.). Tseng and Yun carry these results over to the nonsmooth problem (1) by reformulating it as a smooth problem over the epigraph of the nonsmooth part.

Setting

Throughout, $\mathbb R^n$ carries the Euclidean norm $\|\cdot\|$. The problem is

$$\min_x\; F_c(x) = f(x) + cP(x), \qquad (1)$$

where $c > 0$, $P:\mathbb R^n \to (-\infty,\infty]$ is proper, convex and lower semicontinuous with effective domain $\operatorname{dom}P = \{x \mid P(x) < \infty\}$, and $f$ is continuously differentiable on an open set containing $\operatorname{dom}P$.

The residual at $x \in \operatorname{dom}P$ is

$$d_I(x) = \arg\min_d \Big\{\nabla f(x)^\top d + \tfrac12\|d\|^2 + cP(x+d)\Big\},$$

the unique minimizer of a strongly convex function. A point $x \in \operatorname{dom}P$ is stationary if the one-sided directional derivative satisfies $F_c'(x;d) \ge 0$ for every direction $d$; $\bar X$ denotes the set of stationary points and $\operatorname{dist}(x,\bar X)$ the distance to it. A point is stationary exactly when its residual vanishes.

Assumption 2(a) asks that $\bar X \neq \emptyset$ and that for every level $\zeta$ there be constants $\tau, \epsilon > 0$ with

$$\operatorname{dist}(x,\bar X) \le \tau\|d_I(x)\| \quad\text{whenever } F_c(x) \le \zeta,\ \|d_I(x)\| \le \epsilon.$$

Writing $\operatorname{epi}P = \{(x,\xi) \mid P(x) \le \xi\}$, problem (1) is equivalent to the smooth problem $\min\{f(x) + c\xi \mid (x,\xi) \in \operatorname{epi}P\}$ (41). Its projection residual at $(x,\xi)$ is the optimal solution $(\tilde d,\tilde\delta)$ of

$$\min_{(d,\delta)}\Big\{\nabla f(x)^\top d + \tfrac12\|d\|^2 + \tfrac12\delta^2 + c\delta \;\Big|\; (x+d,\xi+\delta) \in \operatorname{epi}P\Big\}. \qquad (42)$$

$P$ is polyhedral if $\operatorname{epi}P$ is the solution set of finitely many linear inequalities in $(x,\xi)$; $f$ is quadratic if $f(x) = \tfrac12 x^\top Ax + b^\top x + c_0$ with $A$ symmetric, not necessarily positive semidefinite. The four problem classes are:

• C1: $f$ quadratic, $P$ polyhedral;
• C2: $f(x) = g(Ex) + q^\top x$ with $g$ strongly convex and differentiable on $\mathbb R^m$, $\nabla g$ Lipschitz, $P$ polyhedral;
• C3: $f(x) = \max_{y\in Y}\{(Ex)^\top y - g(y)\} + q^\top x$ with $Y$ polyhedral and $g$ as in C2, $P$ polyhedral;
• C4: $f$ strongly convex with $\nabla f$ Lipschitz on $\operatorname{dom}P$ (condition (22)).

Formalization targets

Goal: Theorem 4 (p. 412)

$$\Big(\bar X \neq \emptyset \ \wedge\ (\mathrm{C1} \vee \mathrm{C2} \vee \mathrm{C3})\Big) \ \vee\ \mathrm{C4} \;\Longrightarrow\; \text{Assumption 2(a)}.$$

Under C4, nonemptiness of $\bar X$ is part of the conclusion. The constants $\tau, \epsilon$ are existential, so the goal does not depend on any particular estimate.

Milestone: Lemma 6 (p. 410)

If $P$ is Lipschitz on $\operatorname{dom}P$ with constant $K$, there is $\kappa > 0$ depending only on $K$ with

$$\|(\tilde d,\tilde\delta)\| \le \kappa\,\|d_I(x)\| \qquad (x \in \operatorname{dom}P,\ \xi = P(x)).$$

Milestone: Lemma 7 (p. 411)

If $\bar X \neq \emptyset$ and C1, C2 or C3 holds, then for every $\zeta$ there are $\tau', \epsilon' > 0$ with

$$\operatorname{dist}(x,\bar X) \le \tau'\|(\tilde d,\tilde\delta)\| \quad\text{whenever } F_c(x) \le \zeta,\ \|(\tilde d,\tilde\delta)\| \le \epsilon', \qquad (44)$$

where $(\tilde d,\tilde\delta)$ solves (42) with $\xi = P(x)$.

A further item, not a milestone, states the remark after Assumption 2 (p. 404) that Assumption 2(b) (stationary points with different objective values are uniformly separated) holds whenever $f$ is convex.

Significance

Theorem 4 is what makes the linear convergence theorem of the paper (Theorem 2, the subject of the second mission of this series) applicable. Through C1 it covers every problem with a quadratic $f$ and a polyhedral $P$, including the Lasso and $\ell_1$-regularized or bound-constrained quadratic programs, without convexity of $f$. Through C2 it covers losses of the form $g(Ex)$ with $g$ strongly convex, such as least squares $\tfrac12\|Ex - b\|^2$, where $f$ itself is not strongly convex because $E$ may have a nontrivial kernel. These error bounds became a standard tool for linear rates of proximal and coordinate methods without strong convexity.

The paper's result is proved, but Lemma 7 is proved by citation: it applies three error bounds of Luo and Tseng to the reformulation (41). A complete formalization therefore requires those error bounds for smooth problems over polyhedral sets, which are not in Mathlib. To our knowledge none of the results of this mission has a machine-checked proof.

Difficulty

Lemma 6 and the C4 case of Theorem 4 are short inequality arguments. The difficulty is Lemma 7. The natural first idea, proving an error bound from strong convexity, fails for C1–C3: $f$ may be nonconvex (C1) or have a degenerate Hessian (C2, C3), and the objective $f(x) + c\xi$ of (41) is never strongly convex in $(x,\xi)$. What replaces strong convexity is the polyhedral structure of $\operatorname{epi}P$, and the cited Luo–Tseng error bounds that exploit it are substantial results in their own right, each of which has to be established for a projection residual over a general polyhedron in $\mathbb R^{n+1}$.

Formalization scope

$\mathbb R^n$ is EuclideanSpace ℝ (Fin n), with indices $0,\dots,n-1$. $P$ is the pair $(D, P)$, its effective domain and its finite values there, with "proper, convex, lsc" given by the published ProxNewton.Inexact.IsProperClosedConvex. The value $+\infty$ is never computed: every statement quantifies over $x \in D$, and membership in $\operatorname{epi}P$ is x + d ∈ D ∧ P (x + d) ≤ ξ + δ. $d_I(x)$ is a chosen minimizer of its subproblem, which at $x \in D$ is the paper's unique one. The solutions of (42) enter through a predicate, and every bound is asserted for every optimal solution. Stationarity is the liminf form of $F_c'(x;d) \ge 0$. $\operatorname{dist}$ is the infimum distance. Assumption 2(a) quantifies over every real $\zeta$, which agrees with the paper's "$\zeta \ge \min F_c$" because the condition is vacuous below $\inf F_c$. In Lemma 6 the constant $\kappa$ is quantified before the dimension and the data, so it depends on the Lipschitz constant only.

Explicit choices relative to the printed text:

• In C4, strong convexity and (22) are required on $\operatorname{dom}P$ only, since $f$ is only assumed smooth near $\operatorname{dom}P$. This hypothesis is weaker than the paper's, so the formal theorem is at least as strong.
• In C2 and C3, $E$ is a continuous linear map $\mathbb R^n \to \mathbb R^m$ (equivalently an $m \times n$ matrix). Strong convexity has a positive modulus and $\nabla g$ is globally Lipschitz.
• In C3, the maximum must be attained at every $x$, which forces $Y \neq \emptyset$, as the paper's formula presumes.
• Assumption 2(b) for convex $f$ is stated with convexity on $\operatorname{dom}P$.

A polyhedrality notion that admits only affine or constant $P$, an Assumption 2(a) whose constants are chosen after $x$, and a Lemma 7 whose residual is $d_I(x)$ instead of the solution of (42) would all trivialize or change the statements. The encoding rules out each of them: polyhedral means a finite system of linear inequalities in $(x,\xi)$, the constants precede $x$, and Lemma 7 is stated for (42).

Welcome contributions: proofs of Lemma 6 and of the C4 case; Lipschitz continuity of polyhedral functions on their domain (Rockafellar–Wets, Example 9.35); the Luo–Tseng error bounds for affine variational inequalities and for composite strongly convex objectives over polyhedra, which can be reused well beyond this mission.

Selected references

• P. Tseng, S. Yun, A coordinate gradient descent method for nonsmooth separable minimization, Math. Program. Ser. B 117 (2009) 387–423. https://doi.org/10.1007/s10107-007-0170-0
• Z.-Q. Luo, P. Tseng, Error bounds and the convergence analysis of matrix splitting algorithms for the affine variational inequality problem, SIAM J. Optim. 2 (1992) 43–54. https://doi.org/10.1137/0802004
• Z.-Q. Luo, P. Tseng, On the linear convergence of descent methods for convex essentially smooth minimization, SIAM J. Control Optim. 30 (1992) 408–425. https://doi.org/10.1137/0330025
• Z.-Q. Luo, P. Tseng, On the convergence rate of dual ascent methods for linearly constrained convex minimization, Math. Oper. Res. 18 (1993) 846–867. https://doi.org/10.1287/moor.18.4.846
• Z.-Q. Luo, P. Tseng, Error bounds and convergence analysis of feasible descent methods: a general approach, Ann. Oper. Res. 46 (1993) 157–178. https://doi.org/10.1007/BF02096261
• R. T. Rockafellar, R. J.-B. Wets, Variational Analysis, Springer, 1998. https://doi.org/10.1007/978-3-642-02431-3

Motivation

Many problems in statistics, signal processing and machine learning minimize the sum of a smooth loss and a convex, nonsmooth, separable regularizer: smooth optimization with $\ell_1$-regularization, bound-constrained optimization, the group Lasso and, through duality, support vector regression are special cases listed in §1 of the paper discussed below. For such problems, methods that update one coordinate, or one small block of coordinates, at a time are often the method of choice, because each update is cheap and exploits separability.

Paul Tseng and Sangwoon Yun (Math. Program. Ser. B 117 (2009) 387–423) introduced the coordinate gradient descent (CGD) method for this class: at each iteration it minimizes a quadratic model of the smooth part plus the exact nonsmooth part over a chosen block of coordinates, then takes an Armijo step. Their paper proves global convergence (Theorem 1) and, under a local error bound, linear convergence (Theorems 2 and 3). This mission is the second of a three-mission series on the paper; it targets the linear convergence theorem for the cyclic (Gauss–Seidel-type) choice of blocks.

The analysis follows Luo and Tseng's error-bound framework for smooth constrained problems (Ann. Oper. Res. 46 (1993) 157–178), which established linear rates for feasible descent methods without strong convexity. Tseng and Yun extend that framework to a nonsmooth convex term.

Setting

Let $\Re^n$ carry the Euclidean norm $\|\cdot\|$, let $\mathcal N=\{1,\dots,n\}$, and for $\mathcal J\subseteq\mathcal N$ let $x_{\mathcal J}$ be the subvector of coordinates in $\mathcal J$. The problem (1) is

$$\min_x\ F_c(x)=f(x)+cP(x),$$

where $c>0$, $P:\Re^n\to(-\infty,\infty]$ is proper, convex and lower semicontinuous, and $f$ is continuously differentiable on an open set containing $\operatorname{dom}P$.

Given $x\in\operatorname{dom}P$, a nonempty block $\mathcal J$ and a symmetric positive definite matrix $H$, the direction (6) is

$$d_H(x;\mathcal J)=\arg\min_d\Big\{\nabla f(x)^Td+\tfrac12d^THd+cP(x+d)\ \Big|\ d_j=0\ \forall j\notin\mathcal J\Big\}.$$

The CGD method starts at $x^0\in\operatorname{dom}P$, chooses at each iteration $k$ a block $\mathcal J^k$ and a matrix $H^k$, computes $d^k=d_{H^k}(x^k;\mathcal J^k)$, and sets $x^{k+1}=x^k+\alpha^kd^k$. The Armijo rule takes $\alpha^k$ as the largest element of $\{\alpha^k_{\rm init}\beta^j\}_{j\ge0}$ with

$$F_c(x^k+\alpha^kd^k)\le F_c(x^k)+\alpha^k\sigma\Delta^k,\qquad \Delta^k=\nabla f(x^k)^Td^k+\gamma\,d^{kT}H^kd^k+cP(x^k+d^k)-cP(x^k),$$

with $0<\beta,\sigma<1$ and $0\le\gamma<1$. The restricted Gauss–Seidel rule (12) asks for indices $0=t_0<t_1<\cdots$ (the set $\mathcal T$) such that the blocks used during each cycle $t_i\le k<t_{i+1}$ are pairwise disjoint and cover $\mathcal N$. $P$ is block-separable with respect to $\mathcal J$ if $P(x)=P_{\mathcal J}(x_{\mathcal J})+P_{\mathcal J^C}(x_{\mathcal J^C})$.

The analysis uses three assumptions: $\nabla f$ is $L$-Lipschitz on $\operatorname{dom}P$ (22); Assumption 1, $\bar\lambda I\succeq H^k\succeq\underline\lambda I$ with $0<\underline\lambda\le\bar\lambda$; and Assumption 2. Writing $\bar X$ for the set of stationary points of $F_c$ and $d_I(x)$ for the full direction with $H=I$, Assumption 2 says (a) $\bar X\ne\emptyset$ and the local Lipschitzian error bound $\operatorname{dist}(x,\bar X)\le\tau\|d_I(x)\|$ holds on every level set $\{F_c\le\zeta\}$ wherever $\|d_I(x)\|\le\epsilon$; (b) stationary points with different objective values are at least a fixed distance $\delta>0$ apart.

Formalization targets

Goal: Theorem 2(b)

Under (22), Assumptions 1 and 2, the restricted Gauss–Seidel rule, block-separability of $P$ with respect to every $\mathcal J^k$, and Armijo steps with $\sup_k\alpha^k_{\rm init}\le1$ and $\inf_k\alpha^k_{\rm init}>0$:

$$\{F_c(x^k)\}\downarrow-\infty\quad\text{or}\quad\big(\{F_c(x^k)\}_{\mathcal T}\ \text{Q-linear}\ \text{and}\ \{x^k\}_{\mathcal T}\ \text{R-linear}\big).$$

No rate constant is fixed: the goal asserts the shape of the convergence, so it is not invalidated by sharper constants.

Milestones

Lemma 4 (Hölder dependence of the subproblem's solution on its linear term), Lemma 5(a) (a per-block variational inequality), Lemma 5(b) (an explicit stepsize at which the Armijo test passes), Theorem 1(f) (Armijo stepsizes bounded away from zero, $\Delta^k\to0$, $d^k\to0$), and Theorem 2(a) (the residual $\|d_I(x^k)\|$ at the start of a cycle is bounded by the step lengths of the cycle).

Significance

Theorem 2(b) gives a linear rate for a block coordinate method on a nonsmooth, possibly nonconvex composite problem without strong convexity: the rate follows from the error bound, which holds for a large class of problems studied in §6 of the paper (the third mission of this series).

The theorem is proved in the paper; none of it has been machine-checked, as far as a search of the Prove2Me library shows. Formalizing it requires a reusable treatment of extended-valued convex functions through their effective domains, of the coordinate-restricted proximal subproblem, and of the Armijo rule for composite objectives. The milestones Lemma 4 and Lemma 5 are self-contained facts about this subproblem and are reusable for other proximal coordinate methods.

Difficulty

The classical route to a linear rate, as in Luo and Tseng's smooth analysis, derives $F_c(x^{k+1})-\bar\upsilon\le\tau'\|x^{k+1}-x^k\|^2$ from the error bound, where $\bar\upsilon$ is the limit value. With a nonsmooth $P$ this inequality is not available: the authors say so on p. 405, and work with $-\Delta^k$ in place of the quadratic term. A second obstacle is that the error bound controls the full residual $d_I(x^k)$, while each iteration only computes a block direction at a different point; relating the two over a cycle needs the disjointness of the blocks in the restricted rule and block-separability of $P$. Neither the global analysis nor the error bound alone gives the rate.

Formalization scope

$\Re^n$ is EuclideanSpace ℝ (Fin n), with 0-based coordinates. The extended-valued $P$ is encoded as a pair: its effective domain $D=\operatorname{dom}P$ and its finite values on $D$, with "proper convex lsc" given by the published definition ProxNewton.Inexact.IsProperClosedConvex; $F_c$ is never evaluated off $D$, and every statement carries membership in $D$ explicitly. The direction $d_H(x;\mathcal J)$ is the unique minimizer of (6) when $x\in D$ and $H\succ0$ (chosen by Classical.epsilon), and runs of the method carry the minimizer property itself. The Armijo rule takes the first admissible exponent $j$. $\mathcal T$ is a strictly increasing map $t$ with $t_0=0$. Stationarity is $F_c'(x;d)\ge0$ for all $d$, in liminf form. In Assumption 2, $\operatorname{dist}$ is the infimum distance and "for any $\zeta\ge\min F_c$" is "for every real $\zeta$" (equivalent, since the condition is empty below $\inf F_c$). Q-linear convergence includes convergence of the sequence; R-linear convergence is $\|x^{t_i}-\bar x\|\le Cq^i$ with $q<1$. The rates in the goal are along $\mathcal T$, as the paper states.

Explicit choices relative to the printed text:

• Theorem 2(a) is corrected. As printed, $\|d_I(x^k)\|\le\sup_j\alpha^j\,C\,r^k$ fails for small constant stepsizes and fails without block-separability (counterexamples in the item's statement). The milestone states $\|d_I(x^k)\|\le\max\{1,\sup_j\alpha^j\}\,C\,r^k$ under block-separability of $P$ with respect to every $\mathcal J^k$ (both are what the paper's proof gives and what Theorem 2(b) supplies), with the iterates in $\operatorname{dom}P$, and with $C$ quantified before the problem data, so it depends only on $n,L,\underline\lambda,\bar\lambda$.
• Lemma 5(b)'s range $0\le\alpha\le\min\{1,2\underline\lambda(1-\sigma+\sigma\gamma)/L\}$ is stated as $0\le\alpha\le1$ and $\alpha L\le2\underline\lambda(1-\sigma+\sigma\gamma)$, so that $L=0$ gives $[0,1]$ as on paper rather than $\{0\}$.
• Lemma 5(a) holds for every decomposition of $P$ in (20), and $\bar x$ ranges over $\operatorname{dom}P_{\mathcal J}$ (outside it the left side is $-\infty$).
• "$\lim F_c(x^k)>-\infty$" is "the values $F_c(x^k)$ are bounded below"; "$\{F_c(x^k)\}\downarrow-\infty$" is "nonincreasing and tending to $-\infty$"; suprema and infima of stepsizes are explicit bounds.

A statement that replaces Assumption 2 by strong convexity, claims a rate along the whole sequence rather than along $\mathcal T$, or asserts the Q-linear contraction without convergence to a limit is not this theorem and is out of scope.

Lemma 1, the claim after (10), Lemma 3 and Theorem 1(a), which the proof also uses, are milestones of the first mission of this series and are not restated here. Theorem 3 (the same conclusion along the whole sequence for the Gauss–Southwell-q rule) is not included: its proof is omitted in the paper. Contributions welcome: proofs of the milestones, and general facts about the coordinate-restricted proximal subproblem (existence and uniqueness of $d_H(x;\mathcal J)$, optimality conditions) that the milestones rely on.

Selected references

• P. Tseng and S. Yun, A coordinate gradient descent method for nonsmooth separable minimization, Mathematical Programming Ser. B 117 (2009) 387–423. https://doi.org/10.1007/s10107-007-0170-0
• Z.-Q. Luo and P. Tseng, Error bounds and convergence analysis of feasible descent methods: a general approach, Annals of Operations Research 46 (1993) 157–178. https://doi.org/10.1007/BF02096261
• J. M. Ortega and W. C. Rheinboldt, Iterative Solution of Nonlinear Equations in Several Variables, Academic Press, 1970, Chap. 9 (Q- and R-linear convergence). https://doi.org/10.1137/1.9780898719468
• R. T. Rockafellar and R. J.-B. Wets, Variational Analysis, Springer, 1998. https://doi.org/10.1007/978-3-642-02431-3