# A Primal-Dual Frank-Wolfe Algorithm for Linear Programming††thanks: This
research was supported in part by a Discovery Grant from the Natural Sciences
and Engineering Research Council (NSERC) of Canada.

Matthew Hough Department of Combinatorics & Optimization, University of
Waterloo, 200 University Ave. W., Waterloo, ON, N2L 3G1, Canada,
<EMAIL_ADDRESS>Stephen A. Vavasis Department of Combinatorics &
Optimization, University of Waterloo, 200 University Ave. W., Waterloo, ON,
N2L 3G1, Canada<EMAIL_ADDRESS>
###### Abstract

We present two first-order primal-dual algorithms for solving saddle point
formulations of linear programs, namely FWLP (Frank-Wolfe Linear Programming)
and FWLP-P. The former iteratively applies the Frank-Wolfe algorithm to both
the primal and dual of the saddle point formulation of a standard-form LP. The
latter is a modification of FWLP in which regularizing perturbations are used
in computing the iterates. We show that FWLP-P converges to a primal-dual
solution with error ${\mathcal{O}}(1/\sqrt{k})$ after $k$ iterations, while no
convergence guarantees are provided for FWLP. We also discuss the advantages
of using FWLP and FWLP-P for solving very large LPs. In particular, we argue
that only part of the matrix $A$ is needed at each iteration, in contrast to
other first-order methods.

##  1 Introduction

In recent years, data science applications have given birth to problems of
very large scale. This poses a problem for mature LP solvers that require
solving a system of linear equations at each iteration. First-order methods
(FoMs) for linear programming aim to solve LPs in such a way that their most
expensive operation at each iteration is the product of a matrix and a vector.
Their goal is to provide an alternative to the practitioner over LP algorithms
such as the simplex method or interior point methods for large-scale problems.

The Frank-Wolfe algorithm [7], also referred to as the conditional gradient
algorithm, is a FoM for minimizing a smooth convex objective function over a
compact convex set. A major benefit of the Frank-Wolfe algorithm is that each
iteration requires only the solution of a linear optimization problem over a
convex constraint set, a problem which can be solved efficiently over many
constraint sets used in practice. It is known that the Frank-Wolfe algorithm
converges at a rate of ${\mathcal{O}}(1/k)$ [12].

The focus of this paper is on finding optimal solutions to linear programs in
standard form, that is, solving the following optimization problem: c^Tx Ax =
b x ≥0, where we assume throughout that (1) has an optimal solution. The dual
linear program associated with (1) is b^Ty A^Ty ≤c.

We propose two first-order primal-dual algorithms for simultaneously solving
(1) and (1) inspired by the Frank-Wolfe algorithm but using the non-standard
step-size of $1/(k+1)$ analyzed in [8]. We call our algorithms FWLP and
FWLP-P. FWLP, first introduced in [11], is derived from iteratively applying
the Frank-Wolfe algorithm to the primal and dual problems of a modified
saddle-point formulation of (1):

$\min_{{\bm{x}}\in\Delta}\max_{{\bm{y}}\in\Gamma}\mathcal{L}({\bm{x}},{\bm{y}}):={\bm{c}}^{T}{\bm{x}}+{\bm{y}}^{T}({\bm{b}}-A{\bm{x}}),$
(1)

where we define

$\Delta=\\{{\bm{x}}\in{\mathbb{R}}^{n}:{\bm{x}}\geq{\bm{0}},{\bm{e}}^{T}{\bm{x}}\leq\xi\\}\mbox{
and }\Gamma=[-\eta,\eta]^{m}.$ (2)

Let $({\bm{x}}^{*},{\bm{y}}^{*})$ denote an optimal primal-dual pair of
solutions to (1) and (1). The parameters $\xi,\eta>0$ are assumed to be chosen
large enough to ensure that $2\xi\geq\lVert{\bm{x}}^{*}\rVert_{1}$ and
$2\eta\geq\lVert{\bm{y}}^{*}\rVert_{\infty}$, thus describing redundant
constraints. Despite their redundance, these constraints are necessary for
compactness of the feasible sets corresponding to the primal and dual
subproblems solved by FWLP. FWLP-P is a modification of the original FWLP
algorithm in which regularizing perturbations are introduced when computing
the iterates. FWLP-P is of theoretical importance to FWLP since we are able to
prove the convergence of FWLP-P but no convergence proof is known yet for
FWLP.

With the ubiquity of very large-scale problems, FoMs for linear programming
have seen interest of late [1, 2, 4, 14]. Such methods typically apply
existing ideas and algorithms from continuous optimization to linear
programming. Our work is related particularly to [9], where Gidel et al. use
the Frank-Wolfe algorithm to solve general convex-concave saddle point
problems:

$\min_{{\bm{x}}\in\mathcal{X}}\max_{{\bm{y}}\in\mathcal{Y}}f({\bm{x}},{\bm{y}}),$
(3)

where $f$ is a smooth convex-concave function and
$\mathcal{X}\times\mathcal{Y}$ is a convex compact set. In analyzing the
convergence of their method, the authors introduce a potential function
bounding the distance to optimality and then show that this potential function
decreases to zero at a given rate. We employ this technique in Section 3 to
analyze the convergence of FWLP-P. The results of [9] do not directly apply to
(1) because their strong-convexity assumption is not satisfied. The FWLP
algorithm is also closely related to the Generalized Fictitious Play Algorithm
proposed by Hammond in her 1984 PhD thesis [10]. Applying Hammond’s algorithm
to (1) yields an algorithm very similar to FWLP, the only difference being
that the update of ${\bm{y}}_{k}$ is based on ${\bm{x}}_{k}$ instead of
${\bm{x}}_{k+1}$ like in FWLP. For a more in-depth analysis of the
similarities between FWLP and Hammond’s Generalized Fictitious Play, see [11,
Chapter 6.2].

The algorithms FWLP and FWLP-P are described in Section 2. Iterations of FWLP
are simpler and faster than those of FWLP-P, but we do not have a convergence
proof of FWLP. Our preliminary computational tests (not reported here)
indicate that the two algorithms converge at comparable rates.

Our convergence analysis of FWLP-P, presented in Section 3, first introduces a
potential function $U_{k}$ which we show that for the iterates of FWLP-P has
distance ${\mathcal{O}}(1/\sqrt{k})$ from zero after $k$ iterations. To
complete the analysis, we show that as $U_{k}\to 0$, FWLP-P and FWLP converge
to a primal-dual solution of (1). The potential $U_{k}$ is similar to a
traditional primal-dual optimality gap as discussed in Section 4. A secondary
contribution of this paper is the discussion in Section 5 of how FWLP and
FWLP-P can be implemented efficiently in a way that only part of the matrix
$A$ is needed at each iteration.

##  2 FWLP and FWLP-P

At each iteration, FWLP performs a Frank-Wolfe update on the primal and the
dual of (1) using the step-size $1/(k+1)$ instead of the standard step-size of
$1/(k+2)$. The step-size $1/(k+1)$ was first analyzed by Freund and Grigas in
[8]. Notably, FWLP uses the information obtained from the primal update,
${\bm{x}}_{k+1}$, in the computation of the dual update:

$\displaystyle{\bm{r}}_{k+1}$
$\displaystyle:={\mathop{\mathrm{argmin}}}_{{\bm{r}}}\left\\{({\bm{c}}-A^{T}{\bm{y}}_{k})^{T}{\bm{r}}:{\bm{r}}\in\Delta\right\\},$
(4) $\displaystyle{\bm{x}}_{k+1}$
$\displaystyle:=\frac{k}{k+1}{\bm{x}}_{k}+\frac{1}{k+1}{\bm{r}}_{k+1},$ (5)
$\displaystyle{\bm{s}}_{k+1}$
$\displaystyle:={\mathop{\mathrm{argmax}}}_{{\bm{s}}}\left\\{({\bm{b}}-A{\bm{x}}_{k+1})^{T}{\bm{s}}:{\bm{s}}\in\Gamma\right\\},$
(6) $\displaystyle{\bm{y}}_{k+1}$
$\displaystyle:=\frac{k}{k+1}{\bm{y}}_{k}+\frac{1}{k+1}{\bm{s}}_{k+1}.$ (7)

It is not hard to see that steps (4) and (6) above can be written in closed
form (for more detail, see [11, Chapter 6]), giving rise to Algorithm 2.1. In
this algorithm and for the remainder of the paper, ${\bm{e}}_{i}$ denotes the
$i$th column of the identity matrix, whose length is determined from the
context.

1:Starting points
${\bm{x}}_{0}\in{\mathbb{R}}^{n}_{+},{\bm{y}}_{0}\in{\mathbb{R}}^{m}$,
constraint data $A$ and $\bm{b}$, and objective ${\bm{c}}$.

2:Parameters: $\xi,\eta>0$ such that $2\lVert{\bm{x}}^{*}\rVert\leq\xi$ and
$2\lVert{\bm{y}}^{*}\rVert\leq\eta$.

3:for $k=1,2,\ldots$ do

4: Determine
$i={\mathop{\mathrm{argmin}}}_{t}[{\bm{c}}-A^{\top}{\bm{y}}_{k}]_{t}$.

5: if $[{\bm{c}}-A^{\top}{\bm{y}}_{k}]_{i}\geq 0$ then

6: Step towards zero in ${\bm{x}}$:
$\displaystyle{\bm{x}}_{k+1}:=\frac{k}{k+1}{\bm{x}}_{k}.$ (8)

7: else

8: Step toward $\xi$ for $x^{(i)}$, otherwise step toward zero for
$x^{(j)},j\neq i$:  $\displaystyle{\bm{x}}_{k+1}$
$\displaystyle=\frac{k}{k+1}{\bm{x}}_{k}+\frac{\xi}{k+1}{\bm{e}}_{i}.$ (9)

9: end if

10: Step towards $\pm\eta$ for each coordinate in ${\bm{y}}$ according to the
sign pattern of ${\bm{b}}-A{\bm{x}}_{k+1}$:  $\displaystyle{\bm{y}}_{k+1}$
$\displaystyle:=\frac{k}{k+1}{\bm{y}}_{k}+\frac{\eta}{k+1}{\mathop{\mathrm{sgn}}}({\bm{b}}-A{\bm{x}}_{k+1}).$
(10)

11:end for

Algorithm 2.1 FWLP: A primal-dual algorithm for (1) based on Frank-Wolfe [11].

FWLP-P adds the regularizing perturbations $\|{\bm{r}}\|^{2}/(2\sqrt{k})$ and
$\|{\bm{s}}\|^{2}/(2\sqrt{k})$ to steps (8)/(9) and (10) in the above
description of FWLP yielding (13) and (15) respectively in Algorithm 2.2. We
analyze this algorithm in the next section.

Note that the computation of ${\bm{s}}_{k+1}$ in FWLP-P is a separable
optimization problem allowing a simple median algorithm for each coordinate
entry of ${\bm{s}}_{k+1}$. The computation of ${\bm{r}}_{k+1}$ is more
elaborate, but it requires only evaluation of the smallest (most negative)
elements of ${\bm{c}}-A^{T}{\bm{y}}_{k}$. We cover this in more detail in
Section 5 but of note here is that FWLP-P retains some of the benefit of FWLP,
namely sub-linear time computation per iteration, since only a part of $A$ is
required. Note also that the formulas (13) and (15) for ${\bm{r}}_{k+1}$ and
${\bm{s}}_{k+1}$ can be written equivalently as projections:

$\displaystyle{\bm{r}}_{k+1}$
$\displaystyle:={\mathrm{proj}}_{\Delta}(\sqrt{k}(A^{T}{\bm{y}}_{k}-{\bm{c}})),$
(11) $\displaystyle{\bm{s}}_{k+1}$
$\displaystyle:={\mathrm{proj}}_{\Gamma}(\sqrt{k}({\bm{b}}-A{\bm{x}}_{k+1})).$
(12)

We discuss how to compute these projections efficiently in Section 5.2.

1:Starting points
${\bm{x}}_{0}\in{\mathbb{R}}^{n}_{+},{\bm{y}}_{0}\in{\mathbb{R}}^{m}$,
constraint data $A$ and $\bm{b}$, and objective ${\bm{c}}$.

2:Parameters: $\xi,\eta>0$ such that $2\lVert{\bm{x}}^{*}\rVert\leq\xi$ and
$2\lVert{\bm{y}}^{*}\rVert\leq\eta$.

3:for $k=1,2,\ldots$ do $\displaystyle{\bm{r}}_{k+1}$
$\displaystyle:={\mathop{\mathrm{argmin}}}_{{\bm{r}}}\left\\{({\bm{c}}-A^{T}{\bm{y}}_{k})^{T}{\bm{r}}+\frac{\|{\bm{r}}\|^{2}}{2\sqrt{k}}:{\bm{r}}\in\Delta\right\\},$
(13) $\displaystyle{\bm{x}}_{k+1}$
$\displaystyle:=\frac{k}{k+1}{\bm{x}}_{k}+\frac{1}{k+1}{\bm{r}}_{k+1},$ (14)
$\displaystyle{\bm{s}}_{k+1}$
$\displaystyle:={\mathop{\mathrm{argmax}}}_{{\bm{s}}}\left\\{({\bm{b}}-A{\bm{x}}_{k+1})^{T}{\bm{s}}-\frac{\|{\bm{s}}\|^{2}}{2\sqrt{k}}:{\bm{s}}\in\Gamma\right\\},$
(15) $\displaystyle{\bm{y}}_{k+1}$
$\displaystyle:=\frac{k}{k+1}{\bm{y}}_{k}+\frac{1}{k+1}{\bm{s}}_{k+1}.$ (16)

4:end for

Algorithm 2.2 FWLP-P: FWLP with perturbations.

##  3 Convergence analysis of FWLP-P

The forthcoming analysis derives a recursion (28) for potential $U_{k}$ below
satisfied by the iterates of FWLP-P. The recursion has two perturbation terms
$\delta_{k+1}$ and $\epsilon_{k+1}$ that are subsequently bounded in (30) and
(32). We show that the quantity $U_{k}$ goes to 0 like $1/\sqrt{k}$.

###  3.1 Deriving the recursion and potential function

Let us introduce

$\delta_{k+1}:={\bm{r}}_{k+2}^{T}({\bm{c}}-A^{T}{\bm{y}}_{k+1})+\frac{1}{2\sqrt{k+1}}\|{\bm{r}}_{k+2}\|^{2}-{\bm{r}}_{k+1}^{T}({\bm{c}}-A^{T}{\bm{y}}_{k+1})-\frac{k}{2(k+1)\sqrt{k}}\|{\bm{r}}_{k+1}\|^{2}.$
(17)

Rearrange this equation:

$\delta_{k+1}-{\bm{r}}_{k+2}^{T}({\bm{c}}-A^{T}{\bm{y}}_{k+1})-\frac{1}{2\sqrt{k+1}}\|{\bm{r}}_{k+2}\|^{2}=-{\bm{r}}_{k+1}^{T}({\bm{c}}-A^{T}{\bm{y}}_{k+1})-\frac{k}{2(k+1)\sqrt{k}}\|{\bm{r}}_{k+1}\|^{2}.$
(18)

Rewrite the first-term on the RHS of (18), that is,
$-{\bm{r}}_{k+1}^{T}({\bm{c}}-A^{T}{\bm{y}}_{k+1})$, using (16):

$\displaystyle-{\bm{r}}_{k+1}^{T}({\bm{c}}-A^{T}{\bm{y}}_{k+1})$
$\displaystyle=-{\bm{r}}_{k+1}^{T}\left({\bm{c}}-A^{T}\left(\frac{k}{k+1}{\bm{y}}_{k}+\frac{1}{k+1}{\bm{s}}_{k+1}\right)\right)$
$\displaystyle=-\frac{k}{k+1}{\bm{r}}_{k+1}^{T}({\bm{c}}-A^{T}{\bm{y}}_{k})-\frac{1}{k+1}{\bm{r}}_{k+1}^{T}({\bm{c}}-A^{T}{\bm{s}}_{k+1}).$
(19)

Substitute (19) for the first term on the RHS in (18), moving the second term
of (19) $-{\bm{r}}_{k+1}^{T}({\bm{c}}-A^{T}{\bm{s}}_{k+1})/(k+1)$ to the LHS,
thereby rewriting (18) as:

$\delta_{k+1}-{\bm{r}}_{k+2}^{T}({\bm{c}}-A^{T}{\bm{y}}_{k+1})-\frac{1}{2\sqrt{k+1}}\|{\bm{r}}_{k+2}\|^{2}+\frac{1}{k+1}{\bm{r}}_{k+1}^{T}({\bm{c}}-A^{T}{\bm{s}}_{k+1})\\\
=-\frac{k}{k+1}{\bm{r}}_{k+1}^{T}({\bm{c}}-A^{T}{\bm{y}}_{k})-\frac{k}{2(k+1)\sqrt{k}}\|{\bm{r}}_{k+1}\|^{2}\\\
=\frac{k}{k+1}\left(-{\bm{r}}_{k+1}^{T}({\bm{c}}-A^{T}{\bm{y}}_{k})-\frac{1}{2\sqrt{k}}\|{\bm{r}}_{k+1}\|^{2}\right).$
(20)

The point of this algebra is that the 2nd and 3rd terms of the LHS correspond
to the RHS advanced from $k$ to $k+1$, thus setting up some of the terms of
the recursion.

Next, introduce for $k\geq 2$

$\epsilon_{k+1}:=\frac{k}{k+1}\left(-{\bm{s}}_{k+1}^{T}({\bm{b}}-A{\bm{x}}_{k})+\frac{\sqrt{k+1}}{2k}\|{\bm{s}}_{k+1}\|^{2}+{\bm{s}}_{k}^{T}({\bm{b}}-A{\bm{x}}_{k})-\frac{1}{2\sqrt{k}}\|{\bm{s}}_{k}\|^{2}\right).$
(21)

Rewrite this equation as

$\epsilon_{k+1}+\frac{k}{k+1}\left({\bm{s}}_{k+1}^{T}({\bm{b}}-A{\bm{x}}_{k})-\frac{\sqrt{k+1}}{2k}\|{\bm{s}}_{k+1}\|^{2}\right)=\frac{k}{k+1}\left({\bm{s}}_{k}^{T}({\bm{b}}-A{\bm{x}}_{k})-\frac{1}{2\sqrt{k}}\|{\bm{s}}_{k}\|^{2}\right).$
(22)

Rearranging (14) yields ${\bm{x}}_{k}=(k+1){\bm{x}}_{k+1}/k-{\bm{r}}_{k+1}/k$.
This means that the factor ${\bm{s}}_{k+1}^{T}({\bm{b}}-A{\bm{x}}_{k})$
appearing in the LHS of (22) may be rewritten

$\displaystyle{\bm{s}}_{k+1}^{T}({\bm{b}}-A{\bm{x}}_{k})$
$\displaystyle={\bm{s}}_{k+1}^{T}\left({\bm{b}}-A\left(\frac{k+1}{k}{\bm{x}}_{k+1}-\frac{1}{k}{\bm{r}}_{k+1}\right)\right)$
$\displaystyle=\frac{k+1}{k}{\bm{s}}_{k+1}^{T}({\bm{b}}-A{\bm{x}}_{k+1})-\frac{1}{k}{\bm{s}}_{k+1}^{T}({\bm{b}}-A{\bm{r}}_{k+1}).$
(23)

Thus, the LHS of (22) is rewritten:

LHS of (22)
$\displaystyle=\epsilon_{k+1}+\frac{k}{k+1}\left(\frac{k+1}{k}{\bm{s}}_{k+1}^{T}({\bm{b}}-A{\bm{x}}_{k+1})-\frac{1}{k}{\bm{s}}_{k+1}^{T}({\bm{b}}-A{\bm{r}}_{k+1})-\frac{\sqrt{k+1}}{2k}\|{\bm{s}}_{k+1}\|^{2}\right)$
$\displaystyle=\epsilon_{k+1}+{\bm{s}}_{k+1}^{T}({\bm{b}}-A{\bm{x}}_{k+1})-\frac{1}{k+1}{\bm{s}}_{k+1}^{T}({\bm{b}}-A{\bm{r}}_{k+1})-\frac{1}{2\sqrt{k+1}}\|{\bm{s}}_{k+1}\|^{2}.$
(24)

Thus, we have rewritten (22) as

$\epsilon_{k+1}+{\bm{s}}_{k+1}^{T}({\bm{b}}-A{\bm{x}}_{k+1})-\frac{1}{k+1}{\bm{s}}_{k+1}^{T}({\bm{b}}-A{\bm{r}}_{k+1})-\frac{1}{2\sqrt{k+1}}\|{\bm{s}}_{k+1}\|^{2}\\\
=\frac{k}{k+1}\left({\bm{s}}_{k}^{T}({\bm{b}}-A{\bm{x}}_{k})-\frac{1}{2\sqrt{k}}\|{\bm{s}}_{k}\|^{2}\right).$
(25)

Here, we see the correspondence between the 2nd and 4th term on the LHS with
the two terms on the RHS (with $k$ advanced by 1).

Multiply (14) by ${\bm{c}}$ and rearrange to obtain

${\bm{c}}^{T}{\bm{x}}_{k+1}-\frac{1}{k+1}{\bm{c}}^{T}{\bm{r}}_{k+1}=\frac{k}{k+1}{\bm{c}}^{T}{\bm{x}}_{k}.$
(26)

Similarly, from (16),

$-{\bm{b}}^{T}{\bm{y}}_{k+1}+\frac{1}{k+1}{\bm{b}}^{T}{\bm{s}}_{k+1}=-\frac{k}{k+1}{\bm{b}}^{T}{\bm{y}}_{k}.$
(27)

Now add (20), (25), (26), and (27), noting that many quantities on the LHS
cancel, while all quantities on the RHS contain the factor $k/(k+1)$, to
obtain a recursion:

$\delta_{k+1}+\epsilon_{k+1}+U_{k+1}=\frac{k}{k+1}U_{k},$ (28)

where for $k\geq 2$

$U_{k}:=-{\bm{r}}_{k+1}^{T}({\bm{c}}-A^{T}{\bm{y}}_{k})-\frac{1}{2\sqrt{k}}\|{\bm{r}}_{k+1}\|^{2}+{\bm{s}}_{k}^{T}({\bm{b}}-A{\bm{x}}_{k})-\frac{1}{2\sqrt{k}}\|{\bm{s}}_{k}\|^{2}+{\bm{c}}^{T}{\bm{x}}_{k}-{\bm{b}}^{T}{\bm{y}}_{k}.$
(29)

We call $U_{k}$ the potential function, and we will use this in our
convergence analysis to bound the distance from optimality of (1) and (1).

###  3.2 Bounding the potential function

Our next task is to lower bound $\delta_{k+1}$ and $\epsilon_{k+1}$ so that we
can use (28) to develop a more useful bound on the potential function.

######  Lemma 1.

Recall $\epsilon_{k+1}$ defined in (21). We have the bound

$\epsilon_{k+1}\geq-\frac{m\eta^{2}}{6k^{2}\sqrt{k-1}}.$

###### Proof.

By adding and subtracting multiples of $\|{\bm{s}}_{k+1}\|^{2}$ and and
$\|{\bm{s}}_{k}\|^{2}$ inside (21), we obtain

$\epsilon_{k+1}=\frac{k}{k+1}(\epsilon_{k+1}^{\prime}+\epsilon_{k+1}^{\prime\prime}+\epsilon_{k+1}^{\prime\prime\prime}),$

where

$\displaystyle\epsilon_{k+1}^{\prime}$
$\displaystyle:=-{\bm{s}}_{k+1}^{T}({\bm{b}}-A{\bm{x}}_{k})+\frac{1}{2\sqrt{k-1}}\|{\bm{s}}_{k+1}\|^{2}+{\bm{s}}_{k}^{T}({\bm{b}}-A{\bm{x}}_{k})-\frac{1}{2\sqrt{k-1}}\|{\bm{s}}_{k}\|^{2},$
$\displaystyle\epsilon_{k+1}^{\prime\prime}$
$\displaystyle:=\left(\frac{1}{2\sqrt{k-1}}-\frac{1}{2\sqrt{k}}\right)\|{\bm{s}}_{k}\|^{2},$
$\displaystyle\epsilon_{k+1}^{\prime\prime\prime}$
$\displaystyle:=\left(\frac{\sqrt{k+1}}{2k}-\frac{1}{2\sqrt{k-1}}\right)\|{\bm{s}}_{k+1}\|^{2}.$

We see that $\epsilon_{k+1}^{\prime}\geq 0$ because ${\bm{s}}_{k}$ is the
maximizer of
${\bm{s}}^{T}({\bm{b}}-A{\bm{x}}_{k})-\|{\bm{s}}\|^{2}/(2\sqrt{k-1})$
according to the definition (15), while ${\bm{s}}_{k+1}$ is some other
feasible point, and $\epsilon_{k+1}^{\prime}$ is the difference between the
two objective values. We also see that $\epsilon_{k+1}^{\prime\prime}\geq 0$.

As for $\epsilon_{k+1}^{\prime\prime\prime}$, use the estimate
$\|{\bm{s}}_{k}\|\leq\sqrt{m}\eta$ according to (15). Finally, by taking a
common denominator and observing that $(k-1/(3k))^{2}\leq k^{2}-1$, we obtain
a lower bound of $-1/(6k^{2}\sqrt{k-1})$ on the parenthesized factor. Thus,
adding the three contributions,

$\epsilon_{k+1}\geq-\frac{m\eta^{2}}{6k^{2}\sqrt{k-1}}\cdot\frac{k}{k+1}\geq-\frac{m\eta^{2}}{6k^{2}\sqrt{k-1}}.$
(30)

∎

######  Lemma 2.

Recall $\delta_{k+1}$ defined in (17). We have the bound

$\delta_{k+1}\geq-D/k^{3/2},$

where $D$ is a positive constant depending on the data defined in (33) below.

###### Proof.

Split $\delta_{k+1}$ into four terms:

$\delta_{k+1}=\delta_{k+1}^{\prime}+\delta_{k+1}^{\prime\prime}+\delta_{k+1}^{\prime\prime\prime}+\delta_{k+1}^{\mathrm{iv}}$

where

$\displaystyle\delta_{k+1}^{\prime}$
$\displaystyle:={\bm{r}}_{k+2}^{T}({\bm{c}}-A^{T}{\bm{y}}_{k})+\frac{1}{2\sqrt{k}}\|{\bm{r}}_{k+2}\|^{2}-{\bm{r}}_{k+1}^{T}({\bm{c}}-A^{T}{\bm{y}}_{k})-\frac{1}{2\sqrt{k}}\|{\bm{r}}_{k+1}\|^{2},$
$\displaystyle\delta_{k+1}^{\prime\prime}$
$\displaystyle:=({\bm{r}}_{k+2}-{\bm{r}}_{k+1})^{T}A^{T}({\bm{y}}_{k}-{\bm{y}}_{k+1}),$
$\displaystyle\delta_{k+1}^{\prime\prime\prime}$
$\displaystyle:=\left(\frac{1}{2\sqrt{k+1}}-\frac{1}{2\sqrt{k}}\right)\|{\bm{r}}_{k+2}\|^{2},$
$\displaystyle\delta_{k+1}^{\mathrm{iv}}$
$\displaystyle:=\left(\frac{1}{2\sqrt{k}}-\frac{k}{2(k+1)\sqrt{k}}\right)\|{\bm{r}}_{k+1}\|^{2}.$

We observe that $\delta_{k+1}^{\prime}\geq 0$ because ${\bm{r}}_{k+1}$ is the
minimizer of
${\bm{r}}^{T}({\bm{c}}-A^{T}{\bm{y}}_{k})+\|{\bm{r}}\|^{2}/(2\sqrt{k})$
according to (13), whereas ${\bm{r}}_{k+2}$ is some other feasible point.

Next, we turn to $\delta_{k+1}^{\prime\prime}$, a product of three factors.
Starting on the first factor,

$\displaystyle\|{\bm{r}}_{k+2}-{\bm{r}}_{k+1}\|$
$\displaystyle=\|{\mathrm{proj}}_{\Delta}(\sqrt{k+1}(A^{T}{\bm{y}}_{k+1}-{\bm{c}}))-{\mathrm{proj}}_{\Delta}(\sqrt{k}(A^{T}{\bm{y}}_{k}-{\bm{c}})\|$
$\displaystyle\leq\|\sqrt{k+1}(A^{T}{\bm{y}}_{k+1}-{\bm{c}}))-\sqrt{k}(A^{T}{\bm{y}}_{k}-{\bm{c}})\|$
(31)
$\displaystyle=\|(\sqrt{k+1}-\sqrt{k})(A^{T}{\bm{y}}_{k+1}-{\bm{c}})+\sqrt{k}(A^{T}{\bm{y}}_{k+1}-{\bm{c}}-(A^{T}{\bm{y}}_{k}-{\bm{c}}))\|$
$\displaystyle=\|(\sqrt{k+1}-\sqrt{k})(A^{T}{\bm{y}}_{k+1}-{\bm{c}})+\sqrt{k}A^{T}({\bm{y}}_{k+1}-{\bm{y}}_{k})\|$
$\displaystyle\leq|\sqrt{k+1}-\sqrt{k}|\cdot\|A^{T}{\bm{y}}_{k+1}-{\bm{c}}\|+\sqrt{k}\|A\|\cdot\|{\bm{y}}_{k+1}-{\bm{y}}_{k}\|.$

Here, the first line follows from (11), the second (that is, (31)) from the
fact that the Lipschitz constant of ${\mathrm{proj}}_{C}(\cdot)$ is 1 for any
closed nonempty convex set $C$, the third line adds and subtracts the same
term, and the last line applies the triangle inequality and
submultiplicativity.

Now we use the facts that $\sqrt{k+1}-\sqrt{k}\leq 1/\sqrt{k}$, and, from
(16),

${\bm{y}}_{k+1}-{\bm{y}}_{k}=\frac{1}{k+1}({\bm{s}}_{k+1}-{\bm{y}}_{k}),$

and finally, the bounds $\|{\bm{y}}_{k}\|\leq\sqrt{m}\eta$,
$\|{\bm{s}}_{k}\|\leq\sqrt{m}\eta$ to conclude that
$\|{\bm{y}}_{k+1}-{\bm{y}}_{k}\|\leq 2\sqrt{m}\eta/(k+1)$ and thus

$\|{\bm{r}}_{k+2}-{\bm{r}}_{k+1}\|\leq(1/\sqrt{k})\cdot(\|A\|\sqrt{m}\eta+\|{\bm{c}}\|)+(2/\sqrt{k})\|A\|\sqrt{m}\eta.$

This takes care of the first factor in $\delta_{k+1}^{\prime\prime}$. The
middle factor is bound by $\|A\|$, and the third factor
$\|{\bm{y}}_{k+1}-{\bm{y}}_{k}\|$ is bounded by (see the previous paragraph)
$2\sqrt{m}\eta/(k+1)$. Thus, overall, we obtain

$\displaystyle\delta_{k+1}^{\prime\prime}$
$\displaystyle\geq-\|{\bm{r}}_{k+2}-{\bm{r}}_{k+1}\|\cdot\|A\|\cdot\|{\bm{y}}_{k}-{\bm{y}}_{k+1}\|$
$\displaystyle\geq-C/k^{3/2},$

where $C$ depends on the problem data:

$C\leq 2\|A\|\sqrt{m}\eta(3\|A\|\sqrt{m}\eta+\|{\bm{c}}\|).$

For $\delta_{k+1}^{\prime\prime\prime}$, by finding a common denominator and
then multiplying the resulting fraction by $\sqrt{k}+\sqrt{k+1}$, we obtain

$\delta_{k+1}^{\prime\prime\prime}\geq-1/(4k^{3/2})\cdot\|{\bm{r}}_{k+2}\|^{2}\geq-1/(4k^{3/2})\cdot\xi^{2}.$

Finally, one sees that $\delta_{k+1}^{\mathrm{iv}}\geq 0$. Putting all of
these terms together yields:

$\delta_{k+1}\geq-D/k^{3/2},$ (32)

where

$D:=2\|A\|\sqrt{m}\eta(3\|A\|\sqrt{m}\eta+\|{\bm{c}}\|)+\frac{\xi^{2}}{4}.$
(33)

∎

Combining (28), Lemma 1, and Lemma 2, we now have for all $k\geq 2$

$U_{k+1}\leq\frac{k}{k+1}U_{k}+\frac{D}{k^{3/2}}+\frac{m\eta^{2}}{6k^{2}\sqrt{k-1}}.$

In fact, since $\sqrt{k}\geq 1/\sqrt{k-1}$ for all $k\geq 2$, we can enlarge
$D$ to obtain the bound

$U_{k+1}\leq\frac{k}{k+1}U_{k}+\frac{\bar{D}}{k^{3/2}},$ (34)

where

$\bar{D}:=D+\frac{m\eta^{2}}{6}.$

With this bound, we prove the following bound on $U_{k}$.

######  Theorem 1.

Recall $U_{k}$ defined in (29). We have the bound

$U_{k+1}\leq\frac{F}{\sqrt{k}},$

for all $k\geq 2$, where $F:=\max\\{\sqrt{2}U_{2},6\bar{D}\\}$.

###### Proof.

The $k=2$ case follows immediately from the definition of $F$. Suppose
inductively that $U_{k}\leq F/\sqrt{k}$. Then we check:

$\displaystyle U_{k+1}$ $\displaystyle\leq\frac{k}{k+1}U_{k}+\bar{D}/k^{3/2}$
$\displaystyle\leq\frac{kF}{(k+1)\sqrt{k}}+\bar{D}/k^{3/2}$
$\displaystyle=\frac{k^{2}F+\bar{D}(k+1)}{(k+1)k^{3/2}}$
$\displaystyle=\frac{1}{\sqrt{k+1}}\cdot\frac{k^{2}F+\bar{D}(k+1)}{\sqrt{k+1}\cdot
k^{3/2}}$
$\displaystyle\leq\frac{1}{\sqrt{k+1}}\cdot\frac{k^{2}F+\bar{D}(k+1)}{(\sqrt{k}+1/(3\sqrt{k}))k^{3/2}}$
$\displaystyle=\frac{1}{\sqrt{k+1}}\cdot\frac{k^{2}F+\bar{D}(k+1)}{k^{2}+k/3}$
$\displaystyle\leq\frac{1}{\sqrt{k+1}}\cdot\frac{k^{2}F+kF/3}{k^{2}+k/3}$
$\displaystyle=\frac{1}{\sqrt{k+1}}\cdot F.$

Here, we used (34) for the first line, the induction hypothesis for the second
line, the inequality $\sqrt{k+1}\geq\sqrt{k}+1/(3\sqrt{k})$ for $k\geq 1$ on
the 5th line (easy to confirm by squaring) and the assumption $\bar{D}\leq
F/6$ on the 7th line, which implies for all $k\geq 1$ that $\bar{D}(k+1)\leq
Fk/3$.

Thus, $U_{k}\leq F/\sqrt{k}$. ∎

The perturbation terms from (13) and (15) in FWLP-P were essential to this
analysis. Without them, as in FWLP, there is no useful bound on how much the
primal step ${\bm{r}}_{k+2}$ can differ from ${\bm{r}}_{k+1}$ when the index
of the most-violated constraint changes from one iteration to the next. With
the perturbations, we are able to obtain the inequality (31).

###  3.3 Proving convergence of the potential function implies convergence of
the iterates

We have shown that the potential function $U_{k}$ decreases at a rate of
${\mathcal{O}}(1/\sqrt{k})$ as $k$ increases. It remains to show that the
potential function bounds the distance of the iterates from an optimizer. We
first need the following lemmas. Both lemmas have the same flavor: if an
infeasible point for an LP improves on the optimal objective value, then the
amount of objective improvement is bounded in terms of the amount of
infeasibility.

######  Lemma 3.

Let $p^{*}$ denote the optimal value of (1). Assume $p^{*}$ is finite, i.e.,
(1) is feasible and bounded. Let $({\bm{x}}^{*},{\bm{y}}^{*})$ be an arbitrary
optimizing primal-dual pair of (1) and (1). Then for an arbitrary
$\hat{\bm{x}}\geq{\bm{0}}$,

${\bm{c}}^{T}\hat{\bm{x}}\geq
p^{*}-\|{\bm{y}}^{*}\|_{\infty}\|{\bm{b}}-A\hat{\bm{x}}\|_{1}.$ (35)

###### Proof.

Select an arbitrary optimizing pair $({\bm{x}}^{*},{\bm{y}}^{*})$ for (1). Let
${\bm{k}}:={\bm{b}}-A\hat{\bm{x}}$. Consider the LP,

$\begin{array}[]{rl}\min_{{\bm{x}}}&{\bm{c}}^{T}{\bm{x}}\\\
\mbox{s.t.}&A{\bm{x}}={\bm{b}}-{\bm{k}},\\\ &{\bm{x}}\geq{\bm{0}}.\end{array}$
(36)

This LP is clearly feasible since $\hat{\bm{x}}$ satisfies the constraints. It
is also bounded, as argued by contradiction: If (36) were feasible and
unbounded, then there would exist a certificate of unboundedness, which is a
vector ${\bm{w}}$ such that ${\bm{w}}\geq{\bm{0}}$, $A{\bm{w}}={\bm{0}}$, and
${\bm{c}}^{T}{\bm{w}}<0$. However, such a ${\bm{w}}$ would also certify
unboundedness of (1), which we have already assumed to be bounded.

Therefore, the dual of (36), which is,

$\begin{array}[]{rl}\max_{{\bm{y}}}&({\bm{b}}-{\bm{k}})^{T}{\bm{y}}\\\
\mbox{s.t.}&A^{T}{\bm{y}}\leq{\bm{c}},\end{array}$ (37)

has an optimal solution, say $\hat{\bm{y}}$. Since ${\bm{y}}^{*}$ is also
feasible for (37), we have the following chain of inequalities

$\displaystyle p^{*}$ $\displaystyle={\bm{b}}^{T}{\bm{y}}^{*}$ (by strong
duality of (1))
$\displaystyle=({\bm{b}}-{\bm{k}})^{T}{\bm{y}}^{*}+{\bm{k}}^{T}{\bm{y}}^{*}$
$\displaystyle\leq({\bm{b}}-{\bm{k}})^{T}{\bm{y}}^{*}+\|{\bm{k}}\|_{1}\cdot\|{\bm{y}}^{*}\|_{\infty}$
$\displaystyle\leq({\bm{b}}-{\bm{k}})^{T}\hat{\bm{y}}+\|{\bm{k}}\|_{1}\cdot\|{\bm{y}}^{*}\|_{\infty}$
(since $\hat{\bm{y}}$ maximizes (37))
$\displaystyle\leq{\bm{c}}^{T}\hat{\bm{x}}+\|{\bm{k}}\|_{1}\cdot\|{\bm{y}}^{*}\|_{\infty}.$
(by weak duality between (36) and (37))

Recalling ${\bm{k}}={\bm{b}}-A\hat{\bm{x}}$, the final line in this chain
establishes (35). ∎

######  Lemma 4.

Let $d^{*}$ denote the optimal value of (1). Assume $d^{*}$ is finite, i.e.,
(1) is feasible and bounded. Let $({\bm{x}}^{*},{\bm{y}}^{*})$ be an arbitrary
optimizing primal-dual pair of (1) and (1). Then for an arbitrary
$\hat{\bm{y}}$,

${\bm{b}}^{T}\hat{\bm{y}}\leq d^{*}+\|{\bm{x}}^{*}\|_{1}\cdot l,$ (38)

where $l$ measures the infeasibility of $\hat{\bm{y}}$, that is,

$l:=\max(0,\max_{j}\\{{\bm{e}}_{j}^{T}(A^{T}\hat{\bm{y}}-{\bm{c}})\\}).$ (39)

###### Proof.

This proof is analogous to the previous proof. Let
$({\bm{x}}^{*},{\bm{y}}^{*})$ be an arbitrary primal-dual optimizer of (1).
Consider the dual-form LP given by

$\begin{array}[]{rl}\max_{{\bm{y}}}&{\bm{b}}^{T}{\bm{y}}\\\
\mbox{s.t.}&A^{T}{\bm{y}}\leq{\bm{c}}+l{\bm{e}},\end{array}$ (40)

where $l$ is from (39) and ${\bm{e}}$ denotes the vector of all 1’s. This LP
is clearly feasible since $\hat{\bm{y}}$ satisfies the constraint.
Furthermore, it is bounded as argued by contradiction. If (40) were feasible
and unbounded, there would exist a certificate of unboundedness, that is, a
vector ${\bm{z}}$ such that $A^{T}{\bm{z}}\leq{\bm{0}}$ and
${\bm{b}}^{T}{\bm{z}}>0$. However, this certificate would also certify the
unboundedness of (1), which is assumed to be bounded.

Therefore, the dual of (40), which is

$\begin{array}[]{rl}\min_{{\bm{x}}}&({\bm{c}}+l{\bm{e}})^{T}{\bm{x}}\\\
\mbox{s.t.}&A{\bm{x}}={\bm{b}},\\\ &{\bm{x}}\geq{\bm{0}},\end{array}$ (41)

has an optimal solution which we denote $\hat{\bm{x}}$. Since ${\bm{x}}^{*}$
is also feasible for (41), we have the following chain of inequalities:

$\displaystyle d^{*}$ $\displaystyle={\bm{c}}^{T}{\bm{x}}^{*}$ (by strong
duality of (1))
$\displaystyle=({\bm{c}}+l{\bm{e}})^{T}{\bm{x}}^{*}-l{\bm{e}}^{T}{\bm{x}}^{*}$
$\displaystyle\geq({\bm{c}}+l{\bm{e}})^{T}{\bm{x}}^{*}-l\|{\bm{x}}^{*}\|_{1}$
$\displaystyle\geq({\bm{c}}+l{\bm{e}})^{T}\hat{\bm{x}}-l\|{\bm{x}}^{*}\|_{1}$
(by the optimality of $\hat{\bm{x}}$ in (41))
$\displaystyle\geq{\bm{b}}^{T}\hat{\bm{y}}-l\|{\bm{x}}^{*}\|_{1}$
$\displaystyle\mbox{(by weak duality betwen \eqref{eq:dldu} and
\eqref{eq:pldu})}.$

The final line after rearrangement is (38). ∎

We are now ready to prove the penultimate theorem that shows convergence of
FWLP-P. We show that, assuming a primal-dual optimal LP solution exists and
$\xi,\eta$ have been chosen correctly, the $U_{k}$ plus terms that tend to 0
bound the distance from optimality. Note that an analog of Theorem 2 can be
derived for FWLP using a similar proof. This analog is not presented here
since we are not able to show that $U_{k}\to 0$ for FWLP.

######  Theorem 2.

Suppose $k\geq 2$. Then,

$\displaystyle{\bm{x}}_{k}\geq{\bm{0}},$ (42)
$\displaystyle\|{\bm{b}}-A{\bm{x}}_{k}\|_{1}\leq\frac{2U_{k}}{\eta}+\frac{\xi^{2}}{\eta\sqrt{k}}+\frac{\eta}{\sqrt{k-1}},$
(43)
$\displaystyle\max(0,\max_{j}\\{{\bm{e}}_{j}^{T}(A^{T}{\bm{y}}_{k}-{\bm{c}})\\})\leq\frac{2U_{k}}{\xi}+\frac{\xi}{\sqrt{k}}+\frac{\eta^{2}}{\xi\sqrt{k-1}},$
(44) $\displaystyle{\bm{c}}^{T}{\bm{x}}_{k}-{\bm{b}}^{T}{\bm{y}}_{k}\leq
U_{k},$ (45)

provided

$\xi\geq 2\|{\bm{x}}^{*}\|_{1},\qquad\eta\geq 2\|{\bm{y}}^{*}\|_{\infty}.$
(46)

###### Proof.

It is immediate from the algorithm that (42) holds. Next, notice that
${\bm{s}}_{k}$ is a solution to

$\max_{{\bm{s}}\in\Gamma}\left\\{{\bm{s}}^{T}({\bm{b}}-A{\bm{x}}_{k})-\frac{1}{2\sqrt{k-1}}\|{\bm{s}}\|^{2}\right\\},$

so for arbitrary ${\bm{s}}\in\Gamma$,

$\begin{split}{\bm{s}}_{k}^{T}({\bm{b}}-A{\bm{x}}_{k})-\frac{1}{2\sqrt{k}}\|{\bm{s}}_{k}\|^{2}&\geq{\bm{s}}_{k}^{T}({\bm{b}}-A{\bm{x}}_{k})-\frac{1}{2\sqrt{k-1}}\|{\bm{s}}_{k}\|^{2},\\\
&\geq{\bm{s}}^{T}({\bm{b}}-A{\bm{x}}_{k})-\frac{1}{2\sqrt{k-1}}\|{\bm{s}}\|^{2}.\end{split}$
(47)

Setting ${\bm{s}}={\bm{0}}$ gives

${\bm{s}}_{k}^{T}({\bm{b}}-A{\bm{x}}_{k})-\frac{1}{2\sqrt{k}}\|{\bm{s}}_{k}\|^{2}\geq
0.$

Similarly, ${\bm{r}}_{k+1}$ is a solution to

$\min_{{\bm{r}}}\left\\{{\bm{r}}^{T}({\bm{c}}-A^{T}{\bm{y}}_{k})+\frac{1}{2\sqrt{k}}\|{\bm{r}}\|^{2}\right\\},$
(48)

so for arbitrary ${\bm{r}}\in\Delta$,

${\bm{r}}_{k+1}^{T}({\bm{c}}-A^{T}{\bm{y}}_{k})+\frac{1}{2\sqrt{k}}\|{\bm{r}}_{k+1}\|^{2}\leq{\bm{r}}^{T}({\bm{c}}-A^{T}{\bm{y}}_{k})+\frac{1}{2\sqrt{k}}\|{\bm{r}}\|^{2}.$
(49)

Setting ${\bm{r}}={\bm{0}}$ gives

$-{\bm{r}}_{k+1}^{T}({\bm{c}}-A^{T}{\bm{y}}_{k})-\frac{1}{2\sqrt{k}}\|{\bm{r}}_{k+1}\|^{2}\geq
0.$ (50)

The above two results imply with (29) that
$U_{k}\geq{\bm{c}}^{T}{\bm{x}}_{k}-{\bm{b}}^{T}{\bm{y}}_{k}$, thus
establishing (45).

Setting
${\bm{s}}:=\eta\cdot{\mathop{\mathrm{sgn}}}({\bm{b}}-A{\bm{x}}_{k})\in\Gamma$
in (47) gives the bound

${\bm{s}}_{k}^{T}({\bm{b}}-A{\bm{x}}_{k})-\frac{1}{2\sqrt{k}}\lVert{\bm{s}}_{k}\rVert^{2}\geq\eta\lVert{\bm{b}}-A{\bm{x}}_{k}\rVert_{1}-\frac{1}{2\sqrt{k-1}}m\eta^{2}.$
(51)

Similarly, set ${\bm{r}}:=\xi{\bm{e}}_{j}$ where
$j={\mathop{\mathrm{argmax}}}_{j}\\{{\bm{e}}_{j}^{T}(A^{T}{\bm{y}}_{k}-{\bm{c}})\\}$.
Noting ${\bm{r}}\in\Delta$, we can use (49) to obtain the bound

$-{\bm{r}}_{k+1}^{T}({\bm{c}}-A^{T}{\bm{y}}_{k})-\frac{1}{2\sqrt{k}}\|{\bm{r}}_{k+1}\|^{2}\geq\xi\cdot\max_{j}\\{{\bm{e}}_{j}^{T}(A^{T}{\bm{y}}_{k}-{\bm{c}})\\}-\frac{1}{2\sqrt{k}}\xi^{2}.$

In fact, from (50) the left-hand side above is nonnegative. It follows that we
can strengthen the above bound to

$-{\bm{r}}_{k+1}^{T}({\bm{c}}-A^{T}{\bm{y}}_{k})-\frac{1}{2\sqrt{k}}\|{\bm{r}}_{k+1}\|^{2}\geq\xi\cdot
l-\frac{1}{2\sqrt{k}}\xi^{2},$ (52)

where $l=\max(0,\max_{j}\\{{\bm{e}}_{j}^{T}(A^{T}{\bm{y}}_{k}-{\bm{c}})\\})$.
Using Lemma 3, Lemma 4, (29), (51), and (52), we may write

$\displaystyle U_{k}$
$\displaystyle=-{\bm{r}}_{k+1}^{T}({\bm{c}}-A^{T}{\bm{y}}_{k})-\frac{1}{2\sqrt{k}}\|{\bm{r}}_{k+1}\|^{2}+{\bm{s}}_{k}^{T}({\bm{b}}-A{\bm{x}}_{k})-\frac{1}{2\sqrt{k}}\|{\bm{s}}_{k}\|^{2}+{\bm{c}}^{T}{\bm{x}}_{k}-{\bm{b}}^{T}{\bm{y}}_{k},$
$\displaystyle\geq\xi\cdot
l-\frac{1}{2\sqrt{k}}\xi^{2}+\eta\lVert{\bm{b}}-A{\bm{x}}_{k}\rVert_{1}-\frac{1}{2\sqrt{k-1}}m\eta^{2}+{\bm{c}}^{T}{\bm{x}}_{k}-{\bm{b}}^{T}{\bm{y}}_{k},$
$\displaystyle\geq\xi\cdot
l-\frac{1}{2\sqrt{k}}\xi^{2}+\eta\lVert{\bm{b}}-A{\bm{x}}_{k}\rVert_{1}-\frac{1}{2\sqrt{k-1}}m\eta^{2}+p^{*}-\lVert{\bm{y}}^{*}\rVert_{\infty}\lVert{\bm{b}}-A{\bm{x}}_{k}\rVert_{1}-d^{*}-\lVert{\bm{x}}^{*}\rVert_{1}\cdot
l,$ $\displaystyle\geq\xi\cdot
l-\frac{1}{2\sqrt{k}}\xi^{2}+\eta\lVert{\bm{b}}-A{\bm{x}}_{k}\rVert_{1}-\frac{1}{2\sqrt{k-1}}m\eta^{2}+p^{*}-d^{*}-\frac{\eta}{2}\lVert{\bm{b}}-A{\bm{x}}_{k}\rVert_{1}-\frac{\xi}{2}\cdot
l$ $\displaystyle=\frac{\xi}{2}\cdot
l-\frac{1}{2\sqrt{k}}\xi^{2}+\frac{\eta}{2}\lVert{\bm{b}}-A{\bm{x}}_{k}\rVert_{1}-\frac{1}{2\sqrt{k-1}}m\eta^{2},$

where the fourth line used the assumption (46) and the final line used the
fact that $p^{*}=d^{*}$. We are left with the bound

$\frac{\xi}{2}\cdot l+\frac{\eta}{2}\lVert{\bm{b}}-A{\bm{x}}_{k}\rVert_{1}\leq
U_{k}+\frac{\xi^{2}}{2\sqrt{k}}+\frac{m\eta^{2}}{2\sqrt{k-1}},$

where both terms on the left are nonnegative. It follows that each term must
be individually bounded by the right-hand side above. This establishes (43)
and (44). ∎

######  Theorem 3.

The iterates of FWLP-P converge to an $\epsilon$-optimal solution of (1) and
its dual (1) after ${\mathcal{O}}(1/\epsilon^{2})$ iterations.

###### Proof.

This follows immediately by applying the bound from Theorem 1 in Theorem 2. ∎

##  4 Relating $\bm{U_{k}}$ to the standard primal-dual gap

Consider the general saddle-point problem

$\min_{{\bm{x}}\in\mathcal{X}}\max_{{\bm{y}}\in\mathcal{Y}}f({\bm{x}},{\bm{y}}).$

A standard measure of the optimality gap of an iteration
$({\bm{x}}_{k},{\bm{y}}_{k})$ used in the analysis of many primal-dual
algorithms for saddle-point problems (for example, [6, 9, 15]) is the primal-
dual gap:

$\max\left\\{f({\bm{x}}_{k},\hat{{\bm{y}}}_{k})-f(\hat{{\bm{x}}}_{k},{\bm{y}}_{k}):\hat{\bm{x}}_{k}\in\mathcal{X},\hat{\bm{y}}_{k}\in\mathcal{Y}\right\\}.$
(53)

Obviously, the solution to (53) is

$\displaystyle\hat{{\bm{x}}}_{k}={\mathop{\mathrm{argmin}}}\\{f({\bm{x}},{\bm{y}}_{k}):{\bm{x}}\in\mathcal{X}\\},$
$\displaystyle\hat{{\bm{y}}}_{k}={\mathop{\mathrm{argmax}}}\\{f({\bm{x}}_{k},{\bm{y}}):{\bm{y}}\in\mathcal{Y}\\}.$

It is noted in [2] that for the saddle-point problem associated with LP

$\min_{{\bm{x}}\geq{\bm{0}}}\max_{{\bm{y}}\in{\mathbb{R}}^{m}}{\bm{c}}^{T}{\bm{x}}+{\bm{y}}^{T}({\bm{b}}-A{\bm{x}}),$

the primal-dual gap can be infinite, since the feasible set
${\mathbb{R}}^{n}_{+}\times{\mathbb{R}}^{m}$ is unbounded. This is not an
issue for the modified saddle-point formulation (1) used in our analysis,
since the feasible set $\Delta\times\Gamma$ is bounded by virtue of the
redundant constraints.

Recall the definition of $\mathcal{L}({\bm{x}},{\bm{y}})$ from (1). Let

$\mathcal{M}_{k}:=\max\\{\mathcal{L}({\bm{x}}_{k},{\bm{s}})-\mathcal{L}({\bm{r}},{\bm{y}}_{k}):({\bm{r}},{\bm{s}})\in\Delta\times\Gamma\\},$

which is the specialization of (53) to our setting. We argue that
$\mathcal{M}_{k}$ is a perturbation of $U_{k}$ via the following bound:

$\displaystyle|\mathcal{M}_{k}-U_{k}|$
$\displaystyle=\left|\max_{{\bm{s}}\in\Gamma}[({\bm{b}}-A{\bm{x}}_{k})^{T}{\bm{s}}+{\bm{c}}^{T}{\bm{x}}_{k}]-\min_{{\bm{r}}\in\Delta}[({\bm{c}}-A^{T}{\bm{y}}_{k})^{T}{\bm{r}}+{\bm{b}}^{T}{\bm{y}}_{k}]-U_{k}\right|$
$\displaystyle\leq\left|\max_{{\bm{s}}\in\Gamma}\left[({\bm{b}}-A{\bm{x}}_{k})^{T}{\bm{s}}+{\bm{c}}^{T}{\bm{x}}_{k}-\frac{\|{\bm{s}}\|^{2}}{2\sqrt{k-1}}\right]-\min_{{\bm{r}}\in\Delta}\left[({\bm{c}}-A^{T}{\bm{y}}_{k})^{T}{\bm{r}}+{\bm{b}}^{T}{\bm{y}}_{k}+\frac{\|{\bm{r}}\|^{2}}{2\sqrt{k}}\right]-U_{k}\right|$
$\displaystyle\hphantom{\leq}\quad\mbox{}+\max_{({\bm{r}},{\bm{s}})\in\Delta\times\Gamma}\left|\frac{\|{\bm{s}}\|^{2}}{2\sqrt{k-1}}+\frac{\|{\bm{r}}\|^{2}}{2\sqrt{k}}\right|$
$\displaystyle\leq\left|\left[({\bm{b}}-A{\bm{x}}_{k})^{T}{\bm{s}}_{k}+{\bm{c}}^{T}{\bm{x}}_{k}-\frac{\|{\bm{s}}_{k}\|^{2}}{2\sqrt{k-1}}\right]-\left[({\bm{c}}-A^{T}{\bm{y}}_{k})^{T}{\bm{r}}_{k+1}+{\bm{b}}^{T}{\bm{y}}_{k}+\frac{\|{\bm{r}}_{k+1}\|^{2}}{2\sqrt{k}}\right]-U_{k}\right|$
$\displaystyle\hphantom{\leq}\quad\mbox{}+\frac{m\eta^{2}}{2\sqrt{k-1}}+\frac{\xi^{2}}{2\sqrt{k}}$
$\displaystyle\leq\frac{m\eta^{2}+\xi^{2}}{2\sqrt{k-1}}+{\mathcal{O}}(k^{-3/2}).$

Here, the second line adds and subtracts the same terms and then applies the
triangle inequality for $|\cdot|$. The third line uses the definitions of
${\bm{s}}_{k}$ from (15) and ${\bm{r}}_{k}$ from (13). The fourth line uses
(29), noting that the terms all cancel out except for the difference between a
denominator of $2\sqrt{k-1}$ versus a denominator of $2\sqrt{k}$. This small
remainder is written as ${\mathcal{O}}(k^{-3/2})$ on the fourth line.

##  5 Efficient implementation of FWLP and FWLP-P

A major advantage of FWLP and FWLP-P is their low computational cost per
iteration. Naïve implementations of Algorithms 2.1 and 2.2 have iteration cost
bounded by the cost of a full matrix-vector product. We discuss below how this
can be significantly improved with an efficient implementation.

###  5.1 FWLP

Consider, e.g., the iteration $k=1$. Suppose we store $A{\bm{x}}_{1}$. In
Algorithm 2.1 we can perform an extra step to update $A{\bm{x}}_{2}$ by either
computing

$A{\bm{x}}_{2}=\frac{k}{2}A{\bm{x}}_{1},$

or

$A{\bm{x}}_{2}=\frac{k}{2}A{\bm{x}}_{1}+\frac{\xi}{2}A{\bm{e}}_{i},$

where $i$ is the index of the most violated dual constraint, computed as
${\mathop{\mathrm{argmin}}}\\{c_{i}-{\bm{e}}_{i}^{T}A^{T}{\bm{y}}_{1}\\}$.
Ignoring the cost of computing $i$, such an update runs in
${\mathcal{O}}(m+n)$ operations per iteration (cost of updating ${\bm{x}}_{k}$
and $A{\bm{x}}_{k}$). Note that ${\mathcal{O}}(m+n)$ can be reduced to
${\mathcal{O}}(m)$ if we keep track of the product of scaling factors
$k/(k+1)$ of ${\bm{x}}_{k+1}$ in a separate variable. In contrast, the naïve
implementation of the primal update in FWLP, where one instead performs
matrix-vector products, runs in ${\mathcal{O}}(mn)$ operations per iteration.

The dual update for, e.g., $k=1$, is given by

${\bm{y}}_{2}=\frac{1}{2}{\bm{y}}_{1}+\frac{\eta}{2}{\mathop{\mathrm{sgn}}}({\bm{b}}-A{\bm{x}}_{2}),$

which also runs in ${\mathcal{O}}(m)$ since $A{\bm{x}}_{2}$ has already been
computed in the primal step. The naïve implementation again takes
${\mathcal{O}}(mn)$ operations.

To improve the cost of computing the index $i$, first note that in Algorithm
2.1 we only need to compute ${\bm{e}}_{i}^{T}A^{T}{\bm{y}}_{k}$ for indices
$i$ such that there is a possibility that
$i={\mathop{\mathrm{argmin}}}_{j}\\{c_{j}-{\bm{e}}_{j}^{T}A{\bm{y}}_{k}\\}$.
This could be implemented by a data structure to store the indices $[n]$ in
some order so that only the possibly most infeasible indices need to be
considered at each iteration. In more detail, suppose $j_{k}$ indexes the most
violated dual constraint on iteration $k$, i.e.,
$j_{k}={\mathop{\mathrm{argmin}}}_{j}[{\bm{c}}-A^{T}{\bm{y}}_{k}]_{j}$.
Suppose $j\in[n]$ is some other index. Based on the value of
$[{\bm{c}}-A^{T}{\bm{y}}_{k}]_{j}-[{\bm{c}}-A^{T}{\bm{y}}_{k}]_{j_{k}}$ and
prior knowledge of the stepsize (which tends to $0$ with $k$), one knows in
advance that constraint $j$ could not be the most violated constraint prior to
iteration $k^{\prime}$, where $k^{\prime}$ is a computable index satisfying
$k^{\prime}>k$. Then constraint $j$ does not even have to be considered by the
algorithm on all iterations between $k$ and $k^{\prime}$.

Thus, the algorithm maintains some subset of constraints $\mathcal{S}_{k}$ on
iteration $k$ that need to be examined for possibly being the most violated.
The computation of $A^{T}{\bm{y}}_{k}$ in the primal update thus runs in
${\mathcal{O}}(|\mathcal{S}_{k}|\cdot m)$ plus the cost of updating the data
structure to obtain $\mathcal{S}_{k+1}$. Naturally, the estimate
${\mathcal{O}}(|\mathcal{S}_{k}|\cdot m)$ is further reduced if $A$ is sparse.

The total cost of iteration $k$ for this implementation of FWLP is thus

${\mathcal{O}}(n+\lvert\mathcal{S}_{k}\rvert\cdot m),$

plus the cost of updating the proposed data structure and precomputation. For
large problems, one would expect that $\lvert\mathcal{S}_{k}\rvert\ll n$,
allowing significant speedup to be achieved by using the proposed efficient
implementation of FWLP. And, as mentioned earlier, the “$n$” term may be
dropped with a careful implementation for updating ${\bm{x}}_{k}$, and the
“$\lvert\mathcal{S}_{k}\rvert\cdot m$” term is reduced in the presence of
sparsity.

###  5.2 FWLP-P

The iterations of FWLP-P differ from FWLP in that the steps (13) and (15) are
more computationally involved.

To solve (13), we equivalently consider the projection form (11), which
amounts to solving the quadratic programming problem x 12∥w_0 - x∥^2 e^Tx ≤ξ x
≥0, where ${\bm{w}}_{0}=\sqrt{k}(A^{T}{\bm{y}}_{k}-{\bm{c}})$. Expanding,
rescaling, and dropping constant terms gives the equivalent problem x -w^Tx\+
12∥x∥^2 e^Tx ≤1 x ≥0, where ${\bm{w}}={\bm{w}}_{0}/\xi$. The KKT conditions of
(5.2) are as follows

$\displaystyle-{\bm{w}}+{\bm{x}}+\mu{\bm{e}}-{\bm{z}}$
$\displaystyle={\bm{0}},$ (54) $\displaystyle{\bm{e}}^{T}{\bm{x}}$
$\displaystyle\leq 1,$ (55) $\displaystyle\mu({\bm{e}}^{T}{\bm{x}}-1)$
$\displaystyle=0,$ (56) $\displaystyle{\bm{z}}^{T}{\bm{x}}$ $\displaystyle=0,$
(57) $\displaystyle{\bm{x}},\mu,{\bm{z}}$ $\displaystyle\geq{\bm{0}}.$ (58)

We first prove the following lemma, which will become useful once we define
our algorithm for solving (5.2).

######  Lemma 5.

Suppose $({\bm{x}},\mu,{\bm{z}})$ form a KKT solution for (5.2) with
${\bm{e}}^{T}{\bm{x}}=1$. There must exist an index $j\in[n]$ such that
$w_{j}>\mu$.

Furthermore, we have the following cases regardless of whether
${\bm{e}}^{T}{\bm{x}}=1$.

  1. (i)

If $w_{j}>\mu$ for some $j\in[n]$ , then $x_{j}>0$ and $z_{j}=0$.

  2. (ii)

If $w_{j}<\mu$ for some $j\in[n]$, then $x_{j}=0$ and $z_{j}>0$.

  3. (iii)

If $w_{j}=\mu$ for some $j\in[n]$, then $x_{j}=z_{j}=0$.

###### Proof.

By (54), we have

$w_{i}=x_{i}-z_{i}+\mu$ (59)

for all $i\in[n]$. The condition ${\bm{e}}^{T}{\bm{x}}=1$ along with the
nonnegativity of ${\bm{x}}$ from (58) implies that there must exist an index
$j\in[n]$ such that $x_{j}>0$. Now, the complementarity condition (57) along
with the nonnegativity of ${\bm{z}}$ from (58) imply that $z_{j}=0$. It
follows from (59) that $w_{j}=x_{j}+\mu>\mu$.

Now for (i): suppose, $w_{j}>\mu$ for some $j\in[n]$. By (59),
$x_{j}-z_{j}+\mu>\mu$, which implies $x_{j}>z_{j}$. But since $z_{j}\geq 0$ it
follows that $x_{j}>0$. By applying (57) we see that $z_{j}=0$. The proof for
(ii) is analogous.

For (iii) suppose, $w_{j}=\mu$ for some $j\in[n]$. By (59),
$x_{j}-z_{j}+\mu=\mu$, which implies $x_{j}=z_{j}$. By applying (57) we must
have $x_{j}=z_{j}=0$. ∎

We now state Algorithm 5.1 for solving (5.2). Note that the algorithm computes
the KKT multipliers $\mu,{\bm{z}}$ as well as ${\bm{x}}$ in order to
illustrate its correctness, but in the FWLP-P code, the computation of
${\bm{z}}$ is not needed.

1:Linear term coefficient ${\bm{w}}\in{\mathbb{R}}^{n}$.

2:Sort ${\bm{w}}$ and store result in $\bar{{\bm{w}}}$, along with the
permutation function $\sigma:[n]\to[n]$ which maps the indices of
$\bar{{\bm{w}}}$ back to their original positions in ${\bm{w}}$:
$(\bar{{\bm{w}}},\sigma):=\text{sort}({\bm{w}}).$ (60)

3:Compute the cumulative sum of $\bar{{\bm{w}}}$:
$S:=\text{cumsum}(\bar{{\bm{w}}}).$ (61)

4:for $j=1,2,\ldots,n$ do

5: Compute  $\mu:=\frac{S(j)-1}{j}$ (62)

6: if $\bar{w}_{j}\geq\mu$ and either $j=n$, or $\bar{w}_{j+1}\leq\mu$ then

7: Record index $j$ in the variable $J$.

8: break

9: end if

10:end for

11:if $\mu\geq 0$ then

12: Construct a KKT solution such that ${\bm{e}}^{T}{\bm{x}}=1$:
$\begin{split}&{\bm{x}}:={\bm{0}},\\\
&{\bm{x}}(\sigma(1:J)):=\bar{{\bm{w}}}(1:J)-\mu\cdot{\bm{e}}(1:J),\\\
&{\bm{z}}:=\mu{\bm{e}}-{\bm{w}},\\\
&{\bm{z}}(\sigma(1:J)):={\bm{0}}.\end{split}$ (63)

13: return

14:end if

15:Otherwise, construct a KKT solution such that ${\bm{e}}^{T}{\bm{x}}<1$:
$\begin{split}&{\bm{x}}:=\max({\bm{0}},{\bm{w}}),\\\ &\mu:=0,\\\
&{\bm{z}}:={\bm{x}}-{\bm{w}}.\end{split}$ (64)

16:return

Algorithm 5.1 An efficient algorithm for finding a KKT solution of (5.2).

We now prove that Algorithm 5.1 is finite and correct.

######  Theorem 4.

Algorithm 5.1 finds a KKT solution to (5.2) in finite time.

###### Proof.

It is clear that the algorithm runs in linear time except for the sorting.

Consider the case where $({\bm{x}},\mu,{\bm{z}})$ satisfy the KKT conditions
with ${\bm{e}}^{T}{\bm{x}}=1$. Then by Lemma 5 there must exist an index
$j\in[n]$ such that $w_{j}>\mu$, so $J$ must be well-defined on line 6 of the
algorithm. Moreover, this means that for all $j\in[J],\bar{w}_{j}\geq\mu_{j}$,
and by applying Lemma 5 we get that $x_{j}\geq 0$ and $z_{j}=0$ for all
$j\in[J]$. A similar argument tells us that $\bar{w}_{j}<\mu_{j}$ and thus
$x_{j}=0,z_{j}>0$ for $j\in J^{\prime}:=[n]\setminus[J]$. Using that
${\bm{e}}^{T}{\bm{x}}=1$ and that the entries of ${\bm{x}}$ must be zero
outside of $[J]$, as well as $z_{j}=0$ for all $j\in[J]$, we may multiply (54)
by ${\bm{e}}_{J}^{T}$ and rewrite to obtain:

${\bm{e}}_{J}^{T}{\bm{w}}={\bm{e}}_{J}^{T}{\bm{x}}-{\bm{e}}_{J}^{T}{\bm{z}}+J\mu=1+J\mu,$

where ${\bm{e}}_{J}$ is taken to be the vector with ones in entries $[J]$ and
zero everywhere else. Rearranging gives the formula for $\mu$ used in
Algorithm 5.1:

$\mu=\frac{S(J)-1}{J}$

where we note that $S(J)={\bm{e}}_{J}^{T}{\bm{w}}$. It follows from assumption
of this case ${\bm{e}}^{T}{\bm{x}}=1$ that $\mu\geq 0$ in line 10 of the
algorithm. After equations (63), the algorithm has constructed an ${\bm{x}}$
such that

${\bm{x}}_{J}={\bm{w}}_{J}-\mu{\bm{e}}_{J},\quad{\bm{x}}_{J^{{}^{\prime}}}={\bm{0}},$

and a ${\bm{z}}$ such that

${\bm{z}}_{J}={\bm{0}},\quad{\bm{z}}_{J^{\prime}}=\mu{\bm{e}}_{J^{\prime}}-{\bm{w}}_{J^{\prime}}.$

It is now easy to check that the $({\bm{x}},\mu,{\bm{z}})$ constructed by the
algorithm satisfy the KKT conditions.

Now consider the case where $({\bm{x}},\mu,{\bm{z}})$ satisfy the KKT
conditions with ${\bm{e}}^{T}{\bm{x}}<1$. From (56), $\mu=0$ and clearly the
${\bm{x}}$ and ${\bm{z}}$ defined in (64) satisfy the remaining KKT
conditions. ∎

The cost of Algorithm 5.1 is bounded below by the cost of either the matrix-
vector product $A^{T}{\bm{y}}_{k}$ in computing ${\bm{w}}$ or the sort (60),
which take ${\mathcal{O}}(mn)$ and ${\mathcal{O}}(n\log(n))$ operations
respectively. Meanwhile, (14) takes ${\mathcal{O}}(n)$ operations, so overall
the two run in ${\mathcal{O}}(n+mn)$ operations if we assume $m\geq\log(n)$.

Algorithm 5.1 could be sped up by noticing that in order to construct
${\bm{x}}$,we need only knowledge of the most-violated dual constraints, i.e.,
those for which ${\bm{c}}-A^{T}{\bm{y}}_{k}$ is most negative. Recall from
Section 5.2 that the computation of $A^{T}{\bm{y}}_{k}$ would take
${\mathcal{O}}(\lvert\mathcal{S}_{k}\rvert\cdot m)$ operations, plus the cost
of updating the data structure storing the indices $\mathcal{S}_{k}$. The
sorting computation cost would also be decreased since only the largest
entries would need sorting.

We now consider the dual updates (15) and (16). Step (15) can be solved by
considering the projection form
${\mathrm{proj}}_{\Gamma}(\sqrt{k}({\bm{b}}-A{\bm{x}}_{k+1}))$ and noting that
this has a well known closed-form solution [5, Lemma 6.26]:

${\bm{s}}_{k+1}(i)=\max(-\eta,\min(\eta,\sqrt{k}[{\bm{b}}-A{\bm{x}}_{k+1}]_{i})),$

for each $i\in[m]$. The cost of computing $A{\bm{x}}_{k+1}$ is again
${\mathcal{O}}(\lvert\mathcal{S}_{k}\rvert\cdot m)$ because we can scale
$A{\bm{x}}_{k-1}$ and then add the update $A{\bm{r}}_{k}/k$. The number of
nonzero entries in ${\bm{r}}_{k}$ is $|\mathcal{S}_{k+1}|$, which is the
number of entries that partake in the projection.

##  6 Conclusion

We proposed two primal-dual first-order algorithms, namely FWLP and FWLP-P,
for solving linear programming problems and discussed how both algorithms can
be implemented in such a way that significantly improves their efficiency,
especially for large-scale problems. Our convergence analysis of FWLP-P shows
that the algorithm converges to a primal-dual solution with error
${\mathcal{O}}(1/\sqrt{k})$ after $k$ iterations. Despite this, no convergence
proof is known at this time for the simpler and faster algorithm FWLP, and
analysis of this algorithm is a topic for future research.

Another interesting question is how FWLP and FWLP-P can cope with primal or
dual infeasibility. It should be noted that our proof that $U_{k}\rightarrow
0$ (Theorem 1) did not depend on feasibility nor on the correct choice of
$\xi$ and $\eta$, so monitoring $U_{k}$ cannot diagnose these conditions. We
remark that other first-order algorithms have a theory for infeasibility
detection; see, e.g., [3] and [13]. A related question is how the algorithm
can detect whether the two parameters $\xi$ and $\eta$ of FWLP and FWLP-P have
been chosen correctly.

We showed that FWLP-P converges in the sense that the primal and dual
infeasibility measures tend to 0, as does the duality gap. However, we did not
prove convergence of the iterates. In the case that the LP has multiple
optimizers, an open question is whether the algorithms converge to a
particular optimizer.

## References

  * [1] D. Applegate, M. Díaz, O. Hinder, H. Lu, M. Lubin, B. O’Donoghue, and W. Schudy.  Practical large-scale linear programming using primal-dual hybrid gradient.  In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan, editors, Advances in Neural Information Processing Systems, volume 34, pages 20243–20257. Curran Associates, Inc., 2021.
  * [2] D. Applegate, O. Hinder, H. Lu, and M. Lubin.  Faster first-order primal-dual methods for linear programming using restarts and sharpness.  Mathematical Programming, 201(1):133–184, Sep 2023.
  * [3] David Applegate, Mateo Díaz, Haihao Lu, and Miles Lubin.  Infeasibility detection with primal-dual hybrid gradient for large-scale linear programming.  SIAM Journal on Optimization, 34(1):459–484, 2024.
  * [4] K. Basu, A. Ghoting, R. Mazumder, and Y. Pan.  ECLIPSE: An extreme-scale linear program solver for web-applications.  In Hal Daumé III and Aarti Singh, editors, Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 704–714. PMLR, 13–18 Jul 2020\.
  * [5] Amir Beck.  First-Order Methods in Optimization.  Society for Industrial and Applied Mathematics, Philadelphia, PA, 2017\.
  * [6] Antonin Chambolle and Thomas Pock.  A first-order primal-dual algorithm for convex problems with applications to imaging.  Journal of Mathematical Imaging and Vision, 40(1):120–145, May 2011\.
  * [7] M. Frank and P. Wolfe.  An algorithm for quadratic programming.  Naval Research Logistics Quarterly, 3(1-2):95–110, 1956.
  * [8] Robert M. Freund and Paul Grigas.  New analysis and results for the Frank–Wolfe method.  Mathematical Programming, 155(1):199–230, Jan 2016.
  * [9] G. Gidel, T. Jebara, and S. Lacoste-Julien.  Frank-Wolfe algorithms for saddle point problems.  In Proceedings of the 20th International Conference on Artificial Intelligence and Statistics (AISTATS), 2017.
  * [10] J.H. Hammond.  Solving Asymmetric Variational Inequality Problems and Systems of Equations with Generalized Nonlinear Programming Algorithms.  PhD thesis, Massachussetts Institute of Technology, 1984.
  * [11] Matthew Hough.  Solving saddle point formulations of linear programs with Frank-Wolfe.  Master’s thesis, University of Waterloo, 2023.
  * [12] Martin Jaggi.  Revisiting Frank-Wolfe: Projection-free sparse convex optimization.  In Sanjoy Dasgupta and David McAllester, editors, Proceedings of the 30th International Conference on Machine Learning, volume 28 of Proceedings of Machine Learning Research, pages 427–435, Atlanta, Georgia, USA, 17–19 Jun 2013. PMLR.
  * [13] Tao Jiang, Walaa M Moursi, and Stephen A Vavasis.  Range of the displacement operator of PDHG with applications to quadratic and conic programming.  arXiv preprint arXiv:2309.15009, 2023.
  * [14] X. Li, D. Sun, and KC Toh.  An asymptotically superlinearly convergent semismooth Newton augmented Lagrangian method for linear programming.  SIAM Journal on Optimization, 30(3):2410–2440, 2020.
  * [15] Arkadi Nemirovski.  Prox-method with rate of convergence O(1/t) for variational inequalities with Lipschitz continuous monotone operators and smooth convex-concave saddle point problems.  SIAM Journal on Optimization, 15(1):229–251, 2004.