# Efficient Vertex-Oriented Polytopic Projection for Web-scale Applications

Rohan Ramanath, S. Sathiya Keerthi, Yao Pan, Konstantin Salomatin, Kinjal Basu
Work done while at LinkedIn

###### Abstract

We consider applications involving a large set of instances of projecting
points to polytopes. We develop an intuition guided by theoretical and
empirical analysis to show that when these instances follow certain
structures, a large majority of the projections lie on vertices of the
polytopes. To do these projections efficiently we derive a vertex-oriented
incremental algorithm to project a point onto any arbitrary polytope, as well
as give specific algorithms to cater to simplex projection and polytopes where
the unit box is cut by planes. Such settings are especially useful in web-
scale applications such as optimal matching or allocation problems. Several
such problems in internet marketplaces (e-commerce, ride-sharing, food
delivery, professional services, advertising, etc.), can be formulated as
Linear Programs (LP) with such polytope constraints that require a projection
step in the overall optimization process. We show that in the very recent work
(Basu et al. 2020), the polytopic projection is the most expensive step and
our efficient projection algorithms help in gaining massive improvements in
performance.

##  1 Introduction

Euclidean projection onto a polytope is not only a fundamental problem in
computational geometry (He, Leng, and Li 2006; Rambau and Ziegler 1996), but
also has several real-life applications ranging from decoding of low-density
parity-check codes (Wei and Banihashemi 2017) to generating screening rules
for Lasso (Wang et al. 2013). Recent work (Basu et al. 2020; Boyd et al. 2011)
also uses polytope projection as an operation in solving large optimization
problems. Several such large-scale optimization solvers which use polytopic
projection (Basu et al. 2020; Boyd et al. 2011; Parikh and Boyd 2014) are
being used to tackle different problems in the internet industry such as
matching (Azevedo and Weyl 2016; Makari et al. 2013) and multi-objective
optimization (Agarwal et al. 2011, 2012). Other examples include problems in
natural language processing like structured prediction (Smith 2011; Rush and
Collins 2012; Martins et al. 2015), semi-supervised approaches (Chang,
Sundararajan, and Keerthi 2013), multi-class/hierarchical classification
(Keerthi, Sellamanickam, and Shevade 2012), etc.

Given a compact set $\mathcal{C}$ and a point $x$, let $\Pi_{\mathcal{C}}(x)$
denote the (euclidean) projection of $x$ onto $\mathcal{C}$, i.e.,
$\Pi_{\mathcal{C}}(x)=\arg\min_{y\in\mathcal{C}}\|y-x\|$. Certain large scale
optimization problems, especially those arising from applications in the
internet industry, have a huge number of polytope constraint sets,
$\\{\mathcal{C}_{i}\\}$. To solve such problems each iteration of an algorithm
requires the projection of points $\\{x_{i}\\}$ to
$\\{\tilde{x}_{i}=\Pi_{\mathcal{C}_{i}}(x_{i})\\}$. Computing such sets of
projections is usually the key bottleneck (Wei and Banihashemi 2017) and hence
efficiently solving the projections is crucially important to solve such
large-scale optimization problems.

The set of points, $\\{x_{i}\\}$ usually has some structure connecting the
$x_{i}$ that dictates the nature of the projections. Most general purpose
solvers such as ADMM (Boyd et al. 2011), Block Splitting (Parikh and Boyd
2014) and Splitting Conic Solver (O’Donoghue et al. 2019) do not exploit this
special structure. Very recently, Basu et. al. (Basu et al. 2020) proposed
ECLIPSE (shorthand Ecl) - to solve such large scale problems while exploiting
some general overall structure, but still used off-the-shelf projection
solvers such as the one in (Duchi et al. 2008) to tackle the projection step.
As a result, similar to the general purpose algorithms (Boyd et al. 2011;
O’Donoghue et al. 2019; Parikh and Boyd 2014), the major computational
bottleneck for Ecl is also the projection step.

Although there are many different forms of projection structures, in this
paper we focus on a particular structure motivated by applications in the
internet industry. We consider those, for which, a majority of the
projections, $\tilde{x}_{i}$ lies on a vertex of $X_{i}$ and also, over all
$i$, the mean dimension of the minimal face (Grünbaum et al. 2003) of $X_{i}$
that contains $\tilde{x}_{i}$ is very small. We refer to this as the vertex
oriented structure. Even for special polytopes such as the simplex and box-cut
(a hypercube cut by a single plane) general purpose projection algorithms
(Duchi et al. 2008) are not designed to be efficient for the above situation.
The main goal of this paper is to develop special purpose projection
algorithms that (a) are specially efficient for the vertex oriented structure;
and (b) are applicable to special polytopes such as the simplex, box-cut, and
general polytope forms.

To show the efficacy of these special purpose projection algorithms, we switch
out the projection step of Ecl , and run an ablation study over various
internet marketplace problems. We call this new solver DuaLip, a Dual
decomposition-based Linear Programming framework which shows drastic speedups
over the existing state-of-the-art. Moreover, Ecl only tackled simple polytope
constraints such as the unit box or the simplex, which not only reduced the
generality, but also made the solver tough to use right out of the box. By
replacing the projection component with the novel algorithms developed in this
paper, DuaLip is now capable of solving a much wider range real-world
applications, much more efficiently.

The rest of the paper is organized as follows; We first introduce the problem,
the motivating applications, and describe the Ecl algorithm. We then focus on
the efficient projection algorithms in detail and specialize the general
solution to special structured polytopes. We empirically show the impact of
these algorithms on web-scale data, before finally concluding with a
discussion. All formal proofs are given in the full paper (Ramanath et al.
2021).

##  2 Problem Setup

We begin with some notation. Let $x\in{}^{K}$ denote the point that we wish to
project onto a compact polytope $\mathcal{C}$. The projection operator
$\Pi_{\mathcal{C}}(x)$ can be written as
$\Pi_{\mathcal{C}}(x)=\mathop{\mathrm{argmin}}_{y\in\mathcal{C}}\|y-x\|$ where
$\|\cdot\|$ denotes the Euclidean distance. Throughout this paper, we consider
the following set of polytopes:

  1. 1.

Box : ${\mathcal{C}}=\\{x:0\leq x_{k}\leq 1\;\forall k\\}$

  2. 2.

Simplex-Eq : ${\mathcal{C}}=\\{x:x_{k}\geq 0\;\forall k,\;\sum_{k}x_{k}=1\\}$

  3. 3.

Simplex-Iq : ${\mathcal{C}}=\\{x:x_{k}\geq 0\;\forall k,\;\sum_{k}x_{k}\leq
1\\}$

  4. 4.

Box-Cut-Eq : ${\mathcal{C}}=\\{x:0\leq x_{k}\leq 1\;\forall
k,\;\sum_{k}x_{k}=\delta\\}$ ($\delta=$ positive integer, $1<\delta<K$)

  5. 5.

Box-Cut-Iq : ${\mathcal{C}}=\\{x:0\leq x_{k}\leq 1\;\forall
k,\;\sum_{k}x_{k}\leq\delta\\}$ ($\delta=$ positive integer, $1<\delta<K$)

  6. 6.

Parity Polytope: ${\mathcal{C}}=co\,\\{v_{r}\\}$ where the $v_{r}$ are binary
vectors with an even number of 1s and $co$ denotes convex hull.

  7. 7.

General Polytope : ${\mathcal{C}}=co\,\\{v:v\in V\\}$ where $V$ is a finite
set of polytope vertices.

Here we use $E,I$ to denote an equality or inequality sign denoting whether we
are interested in the surface or the closed interior of the polytope. Such
polytopes naturally occur as constraints in different optimization
formulations in the internet industry. We motivate the need of such polytopes
through typical recommender system problems.

We denote users by $i=1,\ldots,I$ and items by $k=1,\ldots,K$. Let $x_{ik}$
denote any association between user $i$ and item $k$, and be the variable of
interest. For example, $x_{ik}$ can be the probability of displaying item $k$
to user $i$. The vectorized version is denoted by $x=(x_{1},\ldots,x_{I})$
where $x_{i}=\\{x_{ik}\\}_{k=1}^{K}$. Throughout this paper we consider
problems of the form:

$\displaystyle\min_{x}~{}c^{T}x~{}~{}~{}\text{s.t.}~{}~{}~{}Ax\leq
b,~{}~{}x_{i}\in\mathcal{C}_{i},~{}i\in[I],$ (1)

where, $A_{m\times n}$ is the constraint matrix, $b_{m\times 1}$ and
$c_{n\times 1}$ are the constraint and objective vectors respectively, and
${\mathcal{C}}_{i}$ are compact polytopes. $x\in{}^{n}$ is the vector of
optimization variables, where $n=IK$ and $[I]=\\{1,\ldots,I\\}$.

### Applications

Application |  Objective ($c_{ik}$) |  Constraint ($a_{ik}$) | Projection | $\bm{K}$ | $\bm{n}$
---|---|---|---|---|---
Email Optimization (Basu et al. 2020) |  sessions |  unsubscribes | Box | $100$ | $10b$
Diversity in Network Growth |  connection |  invitation | Box-Cut-Eq | $20k$ | $2t$
Item Matching (Jobs Recommendation) |  job applies |  budget | Box-Cut-Iq | $1m$ | $100t$
Multiple Slot Assignment (Feed) |  engagement |  revenue | Assignment | $10b$ | $1q$
Slate Optimization (Ads Ranking) |  revenue |  budget | Slate | $10b$ | $1q$
Table 1:  A list of applications with constraints, projections and size. Here
we use -Eq,-Iq to denote an equality or inequality sign in the constraints
(see Section 3). Finally the short hand, $k,m,b,t,q$ denotes thousand,
million, billion, trillion and quintillion ($10^{18}$), respectively.

Basu et. al. (Basu et al. 2020) described two major classes of problems,
volume optimization (focused on email/notifications) and optimal matching (for
network growth). We cover a much larger class of problems using a variety of
different polytopic constraints.

Diversity in Network Growth: “Rich getting richer” is a common phenomenon in
building networks of users (Fleder and Hosanagar 2009). Frequent or power
users tend to have a large network and get the most attention, while
infrequent members tend to lose out. To prevent these from happening, and to
improve the diversity in the recommended users, we can frame the problem as a
LP:

$\displaystyle\max_{x}\quad$
$\displaystyle\sum\nolimits_{ik}x_{ik}c_{ik}\qquad\qquad\;\;\;\;\;\text{(Total
Utility)}$ s.t. $\displaystyle\sum\nolimits_{i}x_{ik}a_{ik}\geq
b_{k}~{}\forall k\qquad\text{(Infrequent User Bound)}$
$\displaystyle\sum\nolimits_{k}x_{ik}=\delta_{i},\;\;0\leq x_{ik}\leq 1$

where $c,a$ are utility and invitation models, $b_{k}$ denotes the minimum
number of invitations to the $k$-th infrequent user, and $\delta_{i}$ denotes
the number of recommendations to the $i$-th user. This makes $\mathcal{C}_{i}$
the Box-Cut-Eq polytope.

Item Matching in Marketplace Setting: In many two-sided marketplaces, there
are creators and consumers. Each item created has an associated budget and the
problem is to maximize the utility under budget constraints. For example, in
the ads marketplace, each ad campaign has a budget and we need to distribute
impressions appropriately. For the jobs marketplace, each paid job has a
budget and the impressions need to be appropriately allocated to maximize job
applications. Each of these problems can be written as

$\displaystyle\max_{x}\;\;$ $\displaystyle\sum\nolimits_{ik}x_{ik}c_{ik}\;\;$
s.t. $\displaystyle\sum\nolimits_{i}x_{ik}a_{ik}\leq b_{k}\;\forall k\in[K],$
(3) $\displaystyle\sum\nolimits_{k}x_{ik}\leq\delta_{i}\text{ and }0\leq
x_{ik}\leq 1$

where $c,a$ are models estimating utility and budget revenue, respectively and
$\delta_{i}$ is the maximum number of eligible items to be shown to the $i$-th
member. Here the polytope constraint is Box-Cut-Iq .

Multiple Slot Assignment: This is an extension of the item ranking problem
where the utility of an item depends on the position in which it was shown
(Keerthi and Tomlin 2007). Consider for each request or query $i$, we have to
rank $L_{i}$ items ($\ell=1,\ldots,L_{i}$) in $S_{i}$ slots
($s=1,\ldots,S_{i}$). We want to maximize the associated utility subject to
constraints on the items. Mathematically this bipartite structure of multiple
items and slots can be framed as:

$\displaystyle\max_{x}\;\;$ $\displaystyle\sum\nolimits_{i\ell s}x_{i\ell
s}c_{i\ell s}$ s.t. $\displaystyle\sum\nolimits_{i\ell s}x_{i\ell s}a_{i\ell
s}^{(j)}\leq c_{j}~{}~{}~{}~{}~{}\forall j\in[J]$
$\displaystyle\sum\nolimits_{s}x_{i\ell s}=1,\sum\nolimits_{\ell}x_{i\ell
s}=1\text{ and }x_{i\ell s}\geq 0\;\forall i,\ell,s$

where $c$ and $a^{(j)}$ are the associated utility and $j$-th constraint
values. Note that the projection set $\mathcal{C}_{i}$ captures the need to
show each item and the fact that each slot can contain only one item.

This formulation involves a special ${\mathcal{C}}_{i}$ and hence, special
projections. We avoid this by using a revised formulation that introduces an
index $k=1,\ldots,K$ with each $k$ denoting a set of distinct items (instead
of being a single item) in the set of slots. The optimal assignment is then
choosing one such assignment per request. Thus we have,

$\displaystyle\max_{x}$ $\displaystyle\sum\nolimits_{ik}x_{ik}c_{ik}$ s.t.
$\displaystyle\sum\nolimits_{ik}x_{ik}a_{ik}^{(j)}\leq
c_{j}~{}~{}~{}~{}~{}\forall j\in[J]$
$\displaystyle\sum\nolimits_{k}x_{ik}=1\text{ and }x_{ik}\geq 0\;\forall i,k$
(5)

Slate Optimization: This is a variant of the multiple slot assignment problem,
where each of the $L_{i}$ items are ranked and each index $k$ (slate)
corresponds to a ranked selection of distinct items to be placed in the slots.
This can be set up as a path search on a DAG; see (Keerthi and Tomlin 2007)
for more details. Both the Slate and Assignment problems can be written using
the Simplex-Eq polytopic constraint. Table 1 describes these applications with
respect to the associated models, projection, and problem size.

In the subsequent sections, we describe how we develop new projection
algorithms to handle such larger classes of problems. To demonstrate the value
of these novel algorithms, we mainly work with a proprietary Jobs Matching
Dataset ($D$). Here, we are trying to solve the problem of the form in (2),
where instead of the Box-Cut constraint we consider the simplex constraint:
$\sum\nolimits_{k}x_{ik}=1\text{ and }0\leq x_{ik}\leq 1$ and we have $n=100$
trillion. We rely on a real-world dataset only to represent the scale of the
problem, i.e. the number of coupling and non-coupling constraints seen in
practice. The data does not contain any personally identifiable information,
respects GDPR and does not contain any other user information. To promote
reproducible research we open source (with FreeBSD License)111
https://github.com/linkedin/dualip the solver with all the efficient
projection operators discussed in the paper. There, we report the solver’s
performance on the open-source movie-lens dataset (Harper and Konstan 2015),
which contains user rating of movies to formulate an optimization problem
(taking a similar approach to (Makari et al. 2013)). Movie-lens is a public
domain matching problem of a much smaller scale than seen in internet
marketplaces.

### Solving the LP

For sake of completeness, we briefly describe Ecl algorithm here. To solve
problem (1), Ecl introduces a new perturbed problem

$\min_{x}~{}~{}c^{T}x+\frac{\gamma}{2}x^{T}x~{}~{}\text{s.t.}~{}~{}Ax\leq
b,~{}x_{i}\in\mathcal{C}_{i},i\in[I]$ (6)

where $\gamma>0$ helps to make the dual objective function smooth; $\gamma$ is
kept small to keep the solution close to that of the original problem (1). To
make the problem (6) amenable to first order methods, Ecl considers the
Lagrangian dual,

$\displaystyle
g_{\gamma}(\lambda)=\min_{x\in\mathcal{C}}~{}~{}\left\\{c^{T}x+\frac{\gamma}{2}x^{T}x+\lambda^{T}(Ax-b)\right\\},$
(7)

where $\mathcal{C}=\Pi_{i=1}^{I}\mathcal{C}_{i}$. Now, by strong duality (Boyd
and Vandenberghe 2004), the optimum objective $g_{\gamma}^{*}$ of the dual

$\displaystyle g_{\gamma}^{*}:=\max_{\lambda\geq 0}~{}g_{\gamma}(\lambda)$ (8)

is the minimum of (6). Ecl shows that $\lambda\mapsto g_{\gamma}(\lambda)$ is
differentiable and the gradient is Lipschitz continuous. Moreover, by
Danskin’s Theorem (Danskin 2012) the gradient can be explicitly expressed as,
$\nabla g_{\gamma}(\lambda)=Ax_{\gamma}^{*}(\lambda)-b$ where,

$\displaystyle x_{\gamma}^{*}(\lambda)$
$\displaystyle=\mathop{\mathrm{argmin}}_{x\in\mathcal{C}}~{}~{}\left\\{c^{T}x+\frac{\gamma}{2}x^{T}x+\lambda^{T}(Ax-b)\right\\}$
$\displaystyle=\big{\\{}\Pi_{\mathcal{C}_{i}}[-\tfrac{1}{\gamma}(A_{i}^{T}\lambda+c_{i})]\big{\\}}_{i=1}^{I}$
(9)

where $\Pi_{\mathcal{C}_{i}}(\cdot)$ is the Euclidean projection operator onto
$\mathcal{C}_{i}$, and, $A_{i}$, $c_{i}$ are the parts of $A$ and $c$
corresponding to $x_{i}$. Based on this Ecl uses accelerated gradient method
as the main optimizer to solve the problem in (8). For more details, we refer
to (Basu et al. 2020).

### Cost Profiling

On large datasets such as $D$ we find that Ecl takes days to solve. This
motivates us to understand which parts can be improved to get the biggest
overall speedup. We profile each of the parts to analyze the performance
bottlenecks of Ecl . Table 2 shows the complexity of each term required to
compute the gradient along with the sample time spent per iteration of the
solver. It is clear that the optimizer’s internal update steps only take
$~{}3\%$ of the time while the rest ($97\%$) is used to compute the gradient
(columns 2 to 6 in Table 2). The results are consistent on other datasets as
well when the number of executors is tuned appropriately.

Profiling | $A_{i}^{T}\lambda$ | $\hat{x}_{i}(\lambda)$ | $A_{i}\hat{x}_{i}$ | $A\hat{x},c\hat{x},||\hat{x}||$ | $\lambda^{T}(Ax-b)$ | LBFGS
---|---|---|---|---|---|---
Complexity | $\mathcal{O}(k)$ | $\mathcal{O}(k\log k)$ | $\mathcal{O}(\mu_{i})$ | $\mathcal{O}(I/w)+\mathcal{O}(w)$ | $\mathcal{O}(K)$ | $\mathcal{O}(K)$
Sample time | 10% | 74% | 2% | 8% | 3% | 3%
Table 2: Complexity analysis of different components on dataset $D_{1}$ with
$w=100$ nodes. $k$, $\mu_{i}$ as defined later. (a) In gray is the time spent
inside the optimizer, (b) in white is the time taken to compute the gradient.
Notation as defined in (10)

For the rest of paper, we only focus on the most expensive ($74\%$) step in
the gradient computation, that is the projection operation to compute
$x_{\gamma}^{*}(\lambda)$ as in (9).

##  3 Efficient Projection Algorithms

Recall the projection problem to be solved: for each $i$, we want to find

$(x_{\gamma}^{*}(\lambda))_{i}=\Pi_{\mathcal{C}_{i}}[-\tfrac{1}{\gamma}(A_{i}^{T}\lambda+c_{i})]=\arg\min_{x_{i}\in{\mathcal{C}}_{i}}\|x_{i}-\hat{x}_{i}\|^{2}$
(10)

where $\hat{x}_{i}=-\frac{1}{\gamma}(A_{i}^{T}\lambda+c_{i})$. To simplify
notation, we will leave out the ‘$i$’ unless it is necessary. The cost
profiling shows that the projection step forms the bulk of the overall
algorithm cost and hence, improving its speed would lead to overall
algorithmic speedup. If projections form $\phi$ fraction of the total solution
cost, then by making projections extremely efficient one can hope to get up to
a speedup of $1/(1-\phi)$. We now focus on designing efficient algorithms that
are well suited to marketplace problems. Ecl covers Box and Simplex-Eq . The
Box projection solution is very efficient, but in our case the Simplex-Eq
projection algorithm (Duchi et al. 2008) is not. This is primarily because,
for a large fraction of $i$’s, the projections are at a vertex of
$\mathcal{C}_{i}$, which can be identified in a more efficient way.

The Main Intuition: We define a corral (Chakrabarty, Jain, and Kothari 2014)
as the convex hull of a subset of the vertices of ${\mathcal{C}}_{i}$ with
least cardinality that contains the projection of $\hat{x}_{i}$. Thus, if the
projection is a vertex, then the vertex itself is the corral; if the
projection lies in the relative interior of an edge formed by two vertices,
then the corral is that edge, and so on. The dimension of the corral is one
less than the cardinality of its vertex subset. A vertex is a corral of
dimension zero, an edge is a corral of dimension 1, and so on. Clearly there
exists a corral that contains the projection, $x_{i}^{*}$.

Our efficient approach is based on the intuition that, for small $\gamma$, the
mean dimension of the corral over all $i$ is small. (As already mentioned
below (6), the value of $\gamma$ is chosen to be close to zero.) Since this is
the main basis of the paper, we explain the reasoning behind the intuition. We
begin by giving a simple example in 2d. Then we state a formal result that
says what happens at $\gamma=0$. We follow that by taking the case of Simplex-
Eq and demonstrating how, as $\gamma$ becomes larger than zero, the mean
corral dimension only increases gradually with $\gamma$. Further, we give
strong empirical support for the intuition. After that, the rest of the
section and §4 focus on developing efficient projection algorithms based on
the intuition.

Let us begin by illustrating the intuition using an example of Simplex-Eq with
$K=2$; see Figure 1. From the figure it is clear that, for small values of
$\gamma$, the probability that the projection lies on a vertex (corral of
dimension zero) is high.

Figure 1: For a fixed feasible $\lambda$, the set of possible placements of
$\hat{x}$ (see (10) and below it for the definition of $\hat{x}$) that project
to the relative interior of the simplex edge (corral of dimension = 1) is only
the infinite rectangle (shown using dashed lines) perpendicular to the edge.
Unless $\hat{x}$ is exactly perpendicular to the edge (shown in red), when
$\gamma$ becomes small, $\hat{x}$ crosses out of that infinite rectangle,
after which a vertex becomes the projection.

We now state a formal result for the $\gamma=0$ case. The proof of Theorem 1
is given in the full paper (Ramanath et al. 2021).

######  Theorem 1.

Let: $x\in\mathcal{C}$ be expressed as $Dx\leq d$, $Ex=e$, $p=(b,c,d,e)$,
$x^{*}(p)$ denote the optimal solution of (1), and
$x_{0}^{*}(p,\lambda)=\arg\min_{x\in{\mathcal{C}}}(c+A^{T}\lambda)^{T}x$ (also
same as (9) with $\gamma\rightarrow 0$). For any $x$ define
$\mu(x)=\sum_{i}\mu_{i}/I$ to be the mean corral dimension where $\mu_{i}$ is
the corral dimension of $x_{i}$. Let $B_{\epsilon}(z)$ denote the open ball
with center $z$ and radius $\epsilon$. Given
$\epsilon_{b},\epsilon_{c},\epsilon_{d},\epsilon_{e}>0$, if we define
$\mathcal{B}:=\\{(\tilde{b},\tilde{c},\tilde{d},\tilde{e}):\tilde{b}\in
B_{\epsilon_{b}}(b),\tilde{c}\in B_{\epsilon_{c}}(c),\tilde{d}\in
B_{\epsilon_{d}}(d),\tilde{e}\in B_{\epsilon_{e}}(e)\\},$ then

  1. 1.

$x^{*}(p)$ is unique and $\mu(x^{*}(p))\leq m/I$ for almost all222‘almost all’
is same as ‘except for a set of measure zero’. $p\in\mathcal{B}$.

  2. 2.

$x_{0}^{*}(p,\lambda)$ is unique and $\mu(x_{0}^{*}(p,\lambda))=0$, for almost
all $p\in\mathcal{B}$ and $\lambda\in{}^{m}$.

Note that since each $\mu_{i}\in\mathbb{Z}$, it directly follows from Theorem
1 that $\tilde{I}$, the number of $i$’s with vertex solutions ($\mu_{i}=0$) is
at least $I-m$. Moreover, $K$, which denotes the dimension of the $x_{i}$,
does not affect $\mu$ at all (due to the structure of $\mathcal{C}_{i}$ and
the repeating nature in $i$). This is beneficial since $K$ is very large in
Assignment and Slate projections (see Table 1). Comparing results 1 and 2 of
Theorem 1, we can see that non-vertex solutions start appearing as the
solution approaches $\lambda_{0}=\arg\max_{\lambda\geq 0}g_{0}(\lambda)$, and
even there, it is well bounded. In web-scale marketplace applications $m/I$ is
very small, say, smaller than $0.01$. Thus, vertex projections dominate and
they occur in more than 99% of the $i$’s even as we come close to
$\lambda_{0}$.

A limitation of the theorem is that it gives a bound on $\mu$ only at
$\gamma=0$. On the other hand, in our solution process, we use small
$\gamma>0$. By Lemma 2 of (Basu et al. 2020) (also see (Mangasarian and Meyer
1979; Friedlander and Tseng 2008)), the solution $x^{*}$ of (1) remains
unchanged under the $\gamma$ perturbation for sufficiently small positive
$\gamma$ values (Figure 1 gives the rough intuition), and so part 1 of the
theorem applies for such positive $\gamma$.

Even the above formal result is limited. To better explain what happens as we
gradually increase $\gamma$ above zero, let us analyze the Simplex-Eq case.
Let us take any one $i$, but omit the subscript $i$ to avoid clumsiness. Let
$z=-(A^{T}\lambda+c)$ so that $\hat{x}=z/\gamma$. We invoke Algorithm 2 (given
in the next section) to give a precise analysis. Given a $\lambda$ and hence,
given $\hat{x}$, we can use Algorithm 2 to quantify when a vertex (corral of
dimension 0) will be the projection, when an edge (corral of dimension 1) will
contain the projection, etc. In general position (think of random
infinitesimal perturbation of $c$) there will be a positive separation between
the components of $z$. Without loss of generality let us take the components
of $z$ to be ordered: $z_{1}>z_{2}>z_{3}>\ldots$. With this in place, it is
easy to use Algorithm 2 and derive the following: (a) Projection will be a
vertex if $0<\gamma\leq(z_{1}-z_{2})$; (b) Projection will be on an edge if
$(z_{1}-z_{2})<\gamma\leq(z_{1}-z_{3})+(z_{2}-z_{3})$; (c) and so on. Thus
projections start lying on higher dimensional corrals only slowly as $\gamma$
is increased. There is a graded monotonic increase in the optimal corral
dimension with respect to $\gamma$.

It is possible to give similar details for constraint sets other than Simplex-
Eq ; it is just that they get quite clumsy. In addition to what is given
above, the other key point is that, as we move towards optimality, part 1 of
Theorem 1 says that only a few $i$’s (at most $m$ of them at optimality where
$m$ is the number of constraints in $Ax\leq b$) will have non-vertex
projections. Going by the graded effect of positive $\gamma$ on optimal corral
dimensions, we get useful indications of what small values of $\mu$ we can
expect as we close in on optimality with small (not necessarily extremely
small) $\gamma$ (also see the empirical results in Table 3).

Let us use empirical analysis to further understand what happens for larger
$\gamma$ perturbations. Table 3 studies the behavior of $\mu$ along
optimization paths of $\max_{\lambda}g_{\gamma}(\lambda)$ for a range of
$\gamma$ values, on the representative dataset $D$. Even with large $\gamma$
values, note that $\mu$ is usually small (less than $3$). Obtaining a
theoretical bound on $\mu$ as a function of $\gamma$ is an interesting topic
of future academic research. Such theory can be useful in many optimization
formulations that invoke projection as a substep.

Corral metric | $\gamma_{1}=1$ | $\gamma_{2}=0.1$ | $\gamma_{3}=0.01$
---|---|---|---
$\mu$ | 2.55 | 0.12 | 0.03
$\tilde{I}/I$ | 0.4522 | 0.9232 | 0.9676
Table 3: Behavior of the mean corral dimension, $\mu$ along optimization paths
of $g_{\gamma}$ starting from $\lambda=0$ for three values of $\gamma$ s.t.
$\gamma_{1}>\gamma_{2}>\gamma_{3}$.

### Vertex-Oriented Approach

Guided by the theory and the empirical observations that (a) vertex solutions
occur most frequently and (b) the mean corral size is small, we devise the
powerful approach to projection described in Algorithm 1.

Algorithm 1 DuaLip vertex-first projection

1:Let $v^{0}$ be the vertex of ${\mathcal{C}}$ nearest to $\hat{x}$.

2:Check if $v^{0}$ is the optimal solution to (10).

3:If $v^{0}$ is not optimal, include new vertices, $v^{r}$ and search over
corrals of increasing dimension.

As we will see later for Simplex-Eq and Box-Cut-Eq , the algorithms designed
with this approach as the basis are a lot more efficient for our problem
solution than off-the-shelf algorithms that are known to be generally
efficient for finding the projection (Duchi et al. 2008). Also, for the
Assignment and Slate cases, it is important to mention that no efficient
projection algorithms exist in the literature. The algorithms that we develop
for them here are new and are well-suited for our problem solution.

Before we give our algorithms for the specially structured polytopes in 4, let
us discuss a useful general purpose, vertex-oriented projection algorithm for
General Polytope. Many classical algorithms exist for solving the projection-
to-convex-hull problem:

$\min_{x\in{\mathcal{C}}}\|x-\hat{x}\|^{2}\mbox{ where
}{\mathcal{C}}=co\,\\{v:v\in V\\}$ (11)

The algorithms in (Gilbert 1966; Barr and Gilbert 1969; Meyer and Polak 1970;
Mitchell, Dem’yanov, and Malozemov 1974; Gilbert, Johnson, and Keerthi ; Wolfe
1976; Chakrabarty, Jain, and Kothari 2014) are a few of them. DuaLip can be
made to work well with any one of these algorithms. More than the algorithm
choice, the following three key pieces, which are all our contributions, are
absolutely essential for efficiency.333For example, our initial implementation
of using the Wolfe algorithm (Wolfe 1976; Chakrabarty, Jain, and Kothari 2014)
without optimizations (a) through (c) made the projection 100x slower on a
test problem. (a) While these algorithms allow an arbitrary corral for
initialization, our vertex-oriented approach always starts with the nearest-
vertex single-point corral. (b) Even outside of these algorithms we check if
this single-point corral is already optimal for (11) and, if so, these
algorithms are never called. Note that the computations in (a) and (b) are an
integral part of our Algorithm 1 (steps 1 and 2 there). (c) While these
algorithms are set up for the projected point being the origin, doing a
transformation of the convex hull points to meet this requirement would lead
to a large inefficiency. It is important that we keep the original vertices as
they are and carefully modify all operations of the algorithms efficiently.
The summary of this is that, for General Polytope, DuaLip is not wedded
specifically to any particular general-purpose algorithm but rather to the
problem (11), and, for the class of problems and applications covered in the
paper, we do a lot more to make a chosen algorithm to work efficiently than
that algorithm itself.

### Wolfe’s algorithm for General Polytope

Due to its popularity, we delve into the classic algorithm given by Wolfe
(Wolfe 1976; Chakrabarty, Jain, and Kothari 2014) for (10), which also nicely
instantiates the approach in Algorithm 1 for the General Polytope case. We
briefly outline this algorithm here. At a general iteration of the Wolfe
algorithm, there is a corral $C$ formed using a subset of vertices $U$, and
$x$, the projection of $\hat{x}$ to $C$. Using optimality conditions it is
easy to show that $x$ solves (11) iff

$\min\\{(x-\hat{x})^{T}v:v\in V\backslash U\\}\geq(x-\hat{x})^{T}x$ (12)

If (12) is violated then the inclusion of $v=\arg\min\\{(x-\hat{x})^{T}v:v\in
V\backslash U\\}$ to the vertex set of $C$ to form the new corral, $\tilde{C}$
is guaranteed to be such that $\tilde{x}$, the projection on $\tilde{C}$,
satisfies $\|\tilde{x}-\hat{x}\|^{2}<\|x-\hat{x}\|^{2}$. Since it is a descent
algorithm and the number of corrals are finite, the algorithm terminates
finitely. For more details see the full paper (Ramanath et al. 2021).

Wolfe’s algorithm can be initialized with an arbitrary corral and $x$ as the
projection to that corral. However, to be in line with our approach in
Algorithm 1, we initialize $x$ to be the vertex in $V$ that is closest to
$\hat{x}$. This step is also simple and efficient, requiring just a linear
scan over the vertices. It is also useful to note that the immediately
following check of optimality, (12) is also well in line with our approach in
Algorithm 1; the Wolfe algorithm terminates right there if the projection is a
vertex solution, which we know to be the case for most $i$.

##  4 Special Structured Polytope Projections

Wolfe’s algorithm as detailed in Section 3 is highly efficient when the corral
dimension is small. It can be made even more efficient for problems with a
special structure such as Simplex, Box-Cut and Parity Polytope.

Simplex Projection: For Simplex-Eq , ${\mathcal{C}}=co\,\\{v_{k}\\}_{k=1}^{K}$
where $v_{k}$ is the vector whose $k$-th component is equal to $1$ and all the
other components are equal to $0$. Ecl employs the algorithm from (Duchi et
al. 2008, Algorithm 1) for Simplex-Eq , with complexity $O(K\log K)$. The
algorithm also requires all components of $\hat{x}$ to be available; in the
cases of Assignment and Slate, this is very expensive.

The algorithm of (Duchi et al. 2008) is inefficient for our problem since in
most situations, the nearest vertex of ${\mathcal{C}}$, which can be done in
$O(K)$ effort, is expected to be the projection. We can make it more efficient
guided by the approach in Algorithm 1 and specifically, the Wolfe algorithm
outlined earlier. The steps are given in Algorithm 2. Since
$\|v_{j}-\hat{x}\|^{2}=1+\|\hat{x}\|^{2}-2\hat{x}_{j}$, the first step is
equivalent to finding the index, ${\tilde{k}}$ that has the largest $\hat{x}$
component, which forms the initializing corral. The remaining steps correspond
to sequentially including the next best index (Step 3), checking if it does
not lead to any further improvement (Steps 4 and 5), and if so, stopping with
the projection. Step 3 is worth pointing out in more detail. It is based on
(12). Because all vertices of ${\mathcal{C}}$ are orthogonal to each other and
$x$ is a linear combination of the vertices in $U$, $x^{T}v=0\;\forall v\in
V\backslash U$. Thus, $\arg\min\\{(x-\hat{x})^{T}v:v\in V\backslash
U\\}=\arg\max\\{\hat{x}^{T}v:v\in V\backslash U\\}$, leading to Step 3 of
Algorithm 2.

Simplex-Iq can be efficiently solved by first obtaining $x^{\mbox{box}}$, the
box projection. If $\sum_{k}x_{k}^{\mbox{box}}\leq 1$, then $x^{\mbox{box}}$
is also the Simplex-Iq projection; else, the Simplex-Eq projection is the the
Simplex-Iq projection.

In the normal case where all components of $\hat{x}$ are available, the
incremental determination of the next best vertex in step 3 can be done using
the max heapify algorithm (Cormen et al. 2001). If the algorithm stops in $q$
steps, then the overall complexity is $\mathcal{O}(K+q\log K)$. (van den Berg
and Friedlander 2008) gives a similar algorithm especially for Simplex-Eq ,
but without the geometrical intuition and the power of its extension to
structured outputs.

Simplex-Eq Projection for Structured Outputs: Consider any structured output
setting in which, for each $i$ there is a parameter vector, $\xi_{i}$ and a
possibly large set of configurations indexed by $k=1,\ldots,K$ such that: (a)
for each given configuration $k$, $c_{ik}$ and $a_{ik}^{(j)}\;\forall j$ can
be easily computed using $\xi_{i}$; and (b) the incremental determination of
the next best $k$ according to the value,
$\hat{c}_{ik}=c_{ik}+\sum_{j}\lambda_{j}a_{ik}^{(j)}$ is efficient. Using a
Simplex-Eq setup it is clear that Algorithm 2 efficiently solves the
projection problem for this case. Note that the $\hat{c}_{ik}$, $c_{ik}$,
$a_{ik}^{(j)}$ do not need to be computed for all $k$. Instead of $O(K)$,
complexity is usually polynomial in terms of the dimension of $\xi_{i}$.
Assignment and Slate are special cases of such a setup, with each $k$
corresponding to a certain path on a bipartite graph and a DAG, respectively.

Algorithm 2 Modified Duchi et al. algorithm for Simplex-Eq

1:$k_{1}\leftarrow\arg\max_{k\in[K]}\hat{x}_{k}$, and set:
$\mathcal{K}=\\{k_{1}\\}$

2:do

3: $\tilde{k}=\arg\max_{k\in[K]\backslash\mathcal{K}}\hat{x}_{k}$.
$\triangleright$ next best vertex

4:
$\alpha=\hat{x}_{\tilde{k}}-\frac{1}{|\mathcal{K}|+1}\left(\sum_{r\in\mathcal{K}\cup\\{\tilde{k}\\}}\hat{x}_{r}-1\right)$.

5: Set $\mathcal{K}\leftarrow\mathcal{K}\cup\\{\tilde{k}\\}$.

6:while $\alpha>0$

7:$\theta=\frac{1}{\mathcal{K}}\left(\sum_{r\in\mathcal{K}}\hat{x}_{r}\right)$

8:return $x_{r}^{*}=\hat{x}^{r}-\theta\;\forall r\in\mathcal{K}$,
$x_{r}^{*}=0\;\forall r\not\in\mathcal{K}$

Assignment projection: Let us return to the notations given in the definition
of Multiple-slot Assignment in Section 2. Let $\xi_{i}=\\{c_{i\ell
s},\\{a_{i\ell s}^{(j)}\\}_{j}\\}_{\ell,s}$. Each index $k$ corresponds to one
assignment, i.e., a set $\\{(l_{t},s_{t})\\}_{t=1}^{S_{i}}$, with
$c_{ik}=\sum_{t=1}^{S_{i}}c_{i\ell_{t}s_{t}}$ and
$a_{ik}^{(j)}=\sum_{t=1}^{S_{i}}a_{i\ell_{t}s_{t}}^{(j)}$. Thus, we have
$\hat{x}_{ik}=c_{ik}+\sum_{j=1}^{m}\lambda_{j}a_{ik}^{(j)}=\sum_{t=1}^{S_{i}}\phi_{i\ell_{t}s_{t}}$
where $\phi_{i\ell s}=c_{i\ell s}+\sum_{j=1}^{m}\lambda_{j}a_{i\ell s}^{(j)}$.
Thus, to pick ordered elements of $\\{\hat{x}_{k}\\}_{k}$, we just need to
consider an assignment problem with cost defined by $\phi_{i\ell s}$. Step 1
of Algorithm 2 corresponds to picking the best assignment. Step 3 corresponds
to incrementally choosing the next best assignments. Efficient polynomial time
algorithms for implementing these steps are well known (Bourgeois and Lassalle
1971; Kuhn 1955; Murty 1968).

Slate projection: The specialization of Algorithm 2 for this case is similar
to what we described above for Assignment . Because ranking of items needs to
be obeyed and there are edge costs between consecutive items in a slate,
dynamic programming can be used to give an efficient polynomial time Algorithm
(Keerthi and Tomlin 2007) for Step 1. For incrementally finding of the next
best slate (Step 3), the ideas in (Haubold et al. 2018) can be used to give
efficient polynomial time algorithms. It is clear from these specializations
that Algorithm 2 is powerful and can be potentially applied more generally to
problems with other complex structures.

Box-Cut Projection: For Box-Cut-Eq , it is easy to see that
${\mathcal{C}}=co\,\\{v_{r}\\}_{r=1}^{K\choose p}$ where each $v_{r}$ is a
vector with $p$ components equal to $1$ and all the other components equal to
$0$. Without enumeration of the vertices, it is easy to apply Wolfe’s
algorithm. First, for any given direction $\eta\in R^{K}$ let us discuss the
determination of $\max_{r}\eta^{T}v_{r}$. This is equivalent to finding the
top $p$ components of $\eta$ and setting them to $1$ with all remaining
components set to $0$. Using a heap this can be done efficiently in
$\mathcal{O}(K)$ time. Let $\hat{v}(\eta)$ denote the vertex of the Box-Cut
polytope thus found.

The initializing step of the Wolfe algorithm is the finding of the vertex $x$
nearest to $\hat{x}$. Since
$\|v_{r}-\hat{x}\|^{2}=p+\|\hat{x}\|^{2}-2\hat{x}^{T}v_{r}$ and the first two
terms do not depend on $v_{r}$, this step is equivalent to finding
$\hat{v}(\hat{x})$. This nearest vertex is set as the initializing corral. A
general step of the Wolfe algorithm requires the checking of (12). Since $U$
is a corral and $x$ is the projection of $\hat{x}$ to the affine hull of $U$,
$(x-\hat{x})^{T}u=(x-\hat{x})^{T}x\;\forall u\in U$. Thus (12) can also be
written as

$\min\\{(x-\hat{x})^{T}v:v\in V\\}\geq(x-\hat{x})^{T}x.$ (13)

So, all that we need to do is to set $\eta=(\hat{x}-x)$ and obtain
$\hat{v}(\eta)$ to check (13) as $-\eta^{T}\hat{v}(\eta)\geq-\eta^{T}x$ to
stop the Wolfe algorithm or, if that fails, use $\hat{v}(\eta)$ as the next
vertex for proceeding with the Wolfe algorithm. Since, by arguments of
sections 3, the number of steps of the Wolfe algorithm will be just $1$ for
most $i$’s, the cost of doing Box-Cut-Eq projection per $i$ is just
$\mathcal{O}(K)$. This way of doing projection is much more efficient for use
in our problem than more general algorithms with $\mathcal{O}(K\log K)$
complexity given in the literature (Dai and Fletcher 2006).

Box-Cut-Iq can be done efficiently similar to the Simplex-Iq projection by
using the Box and the Box-Cut-Eq projection.

Parity Polytope Projection: ADMM based LP decoding of LDPC codes (Wei and
Banihashemi 2017) requires projections to a special polytope called the Parity
(aka Check) polytope, which is ${\mathcal{C}}=co\,\\{v_{r}\\}$ where the
$v_{r}$ are binary vectors with an even number of 1s. The projection
algorithms in (Lin et al. 2019; Wei and Banihashemi 2017) are specially
designed for the case, but they are general purpose algorithms that use the
inequalities-based representation of $\mathcal{C}_{i}$ and are not designed to
be efficient for the vertex oriented structure.

The specialization of the Wolfe algorithm, that provides such an efficiency,
is similar to that for Box-Cut-Eq . The key computations are: (a) finding the
vertex $v$ nearest to $\hat{x}$; and (b) checking (13), which requires, for
any given any direction $\eta$, $\hat{v}(\eta)=\arg\max_{v\in V}\eta^{T}v$.
Let us start with (b). $\hat{v}(\eta)$ is a binary vector and we need to find
the elements of $\hat{v}(\eta)$ having the 1s. Since we are maximizing
$\eta^{T}v$, take the positive indices, $P=\\{i:\eta_{i}>0\\}$. If $|P|$ is
even, choose the set of elements of $\hat{v}(\eta)$ with 1s as $P$. If $|P|$
is odd, let $i_{1}$ be the index of the element of $P$ with the smallest
$\eta_{i}$, and $i_{2}$ be the index of the element of the complement of $P$
with the largest $\eta_{i}$. If $\eta_{i_{1}}+\eta_{i_{2}}>0$, include $i_{1}$
and $i_{2}$ in $P$; else, remove $i_{1}$ from $P$. This is the optimal way of
choosing $|P|$ to be even. Then choose the set of elements of $\hat{v}(\eta)$
with 1s as $P$. Thus, the determination of $\hat{v}(\eta)$ requires just a
linear scan of $\\{\eta_{i}\\}$.

Let us now come to the finding of the vertex $v$ nearest to $\hat{x}$, i.e.,
$\arg\min_{v}\|\hat{x}-v\|^{2}$. Now,
$\|\hat{x}-v\|^{2}=\|\hat{x}\|^{2}+\|v\|^{2}-2\hat{x}^{T}v$. Since
$\|\hat{x}\|^{2}$ is a constant in the finding of $v$, and, for binary
vectors, $\|v\|^{2}=e^{T}v$ where $e$ is a vector of all 1s, we get
$\arg\min_{v}\|\hat{x}-v\|^{2}=\arg\max_{v}\eta^{T}v$, where
$\eta=\hat{x}-0.5e$. With $\eta$ set this way, we get
$\arg\min_{v}\|\hat{x}-v\|^{2}$ to be $\hat{v}(\eta)$, whose determination was
described in the previous paragraph. Thus, the finding of the vertex nearest
to $\hat{x}$ also requires just a linear scan of $\\{\eta_{i}\\}$. It is worth
trying the Wolfe algorithm as an alternative to existing algorithms to project
points to parity polytopes in the ADMM based LP decoding of LDPC codes (Wei
and Banihashemi 2017).

##  5 Experiments

We demonstrate the value of the algorithms designed in §3 and §4 empirically
using dataset $D$. We solve the problem, using different number of cluster
nodes ($w$) and report the speedup in Table 4 by comparing with the large
scale Ecl solver. Note that, since Ecl only supported the Simplex-Eq
projection via off-the-shelf projection algorithms (Duchi et al. 2008), we
focus our results on the Simplex-Eq use case on their data (Basu et al. 2020)
to showcase an apples to apples comparison of the overall LP solver.

# nodes ($w$) | 55 | 85 | 165 | 800
---|---|---|---|---
$\left\lceil I/w\right\rceil$ (in $100k$) | 15 | 10 | 5 | 1
Speedup | 6x | 6.5x | 7.5x | 8x
Table 4:  Speedup is the ratio of time taken by the DuaLiP projection
algorithm to that of Ecl (based on (Duchi et al. 2008)).

The speedup measures the improvement in time taken to aggregate $I/w$
projection operations. Note that for large values $I/w$, both methods have a
fixed cost of going over a lot of data instances that reduces the relative
speedup of the DuaLip projection. From Table 4, it is clear that even for a
simple projection type such as Simplex-Eq , employing off-the-shelf algorithms
(Duchi et al. 2008) may not be ideal. Basu et. al.(Basu et al. 2020) show that
commercial solvers (SDPT3, CPLEX, CLP, GUROBI) do not scale to the problem-
sizes in web-scale marketplaces and use a distributed version of the Splitting
Conic Solver (SCS (O’Donoghue et al. 2019)) to establish a baseline for scale.
We benchmark DuaLip against Ecl and SCS in Table 5 on real data to show a $7$x
improvement over the state-of-the-art. Ecl and DuaLip exploit the data
sparsity in the constraint matrix ($A$) often present in marketplace data
unlike SCS and other commercial solvers. Additionally, DuaLip benefits from
efficient projection (§3 and §4). We observe similar large efficiency gains on
multiple applications of comparable size using our incremental projection
algorithm. Tables 4 and 5 show that the speed-up depends less on the dataset
and is more a function of the projection and $I/w$.

Problem | Scale $n$ | Time (hr)
---|---|---
DuaLip | Ecl | SCS
Email Optimization (Basu et al. 2020) $nnz(A)=10I\approx 1B$ | $10^{7}$ | 0.11 | 0.8 | 2.0
$10^{8}$ | 0.16 | 1.3 | >24
$10^{9}$ | 0.45 | 4.0 | >24
Item Matching (2) $nnz(A)=100I\approx 10B$ | $10^{10}$ | 0.62 | 4.5 | >24
$10^{11}$ | 1.1 | 7.2 | >24
$10^{12}$ | 2.60 | 11.9 | >24

Table 5: Run time (hrs) for extreme-scale problems on real data. Here,
$nnz(A)$ denotes the number of non-zero entries in $A$ and all runs are on
Spark $2.3$ clusters with up to $800$ processors.

##  6 Discussion

Through empirical evidence and guided by theory, we have developed novel
algorithms for efficient polytope projections under a vertex oriented
structure. Such structures commonly occur in very large-scale industrial
applications. Using these novel algorithms we were able to achieve a drastic
overall speedup compared to the state-of-the-art while solving such web-scale
Linear Programs. This scalability allows us to expand to a much larger class
of problems with complex structured constraints such as Slate and Assignment .
It also opens up the applicability in semi-supervised learning, structured
prediction on text data, optimal marketplace auctions (Edelman, Ostrovsky, and
Schwarz 2007), and rank order of items. All trends reported in this paper hold
across many web marketplace problems (§2). We provide results from one large-
scale internal dataset ($D$) for clarity of impact and open-sourced the solver
on Github444https://github.com/linkedin/dualip with benchmarks on a public-
domain dataset for reproducibility.

## References

  * Agarwal et al. (2011) Agarwal, D.; Chen, B.-C.; Elango, P.; and Wang, X. 2011.  Click Shaping to Optimize Multiple Objectives.  In _KDD_ , 132–140. New York, NY, USA: ACM.  ISBN 978-1-4503-0813-7.
  * Agarwal et al. (2012) Agarwal, D.; Chen, B.-C.; Elango, P.; and Wang, X. 2012.  Personalized Click Shaping Through Lagrangian Duality for Online Recommendation.  In _SIGIR_ , 485–494. New York, NY, USA: ACM.  ISBN 978-1-4503-1472-5.
  * Azevedo and Weyl (2016) Azevedo, E. M.; and Weyl, E. G. 2016.  Matching markets in the digital age.  _Science_ , 352(6289): 1056–1057.
  * Barr and Gilbert (1969) Barr, R.; and Gilbert, E. 1969.  Some efficient algorithms for a class of abstract optimization problems arising on optimal control.  _IEEE Trans. Aut. Contr._ , AC-14: 640–652.
  * Basu et al. (2020) Basu, K.; Ghoting, A.; Mazumder, R.; and Pan, Y. 2020.  ECLIPSE: An Extreme-Scale Linear Program Solver for Web-Applications.  In _ICML 2020_.
  * Bourgeois and Lassalle (1971) Bourgeois, F.; and Lassalle, J.-C. 1971.  An Extension of the Munkres Algorithm for the Assignment Problem to Rectangular Matrices.  14(12): 802–804.
  * Boyd et al. (2011) Boyd, S.; Parikh, N.; Chu, E.; Peleato, B.; and Eckstein, J. 2011.  Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers.  _Foundations and Trends in Machine Learning_ , 3(1): 1–122.
  * Boyd and Vandenberghe (2004) Boyd, S.; and Vandenberghe, L. 2004.  _Convex Optimization_.  Cambridge university press.
  * Chakrabarty, Jain, and Kothari (2014) Chakrabarty, D.; Jain, P.; and Kothari, P. 2014.  Provable Submodular Minimization using Wolfe’s Algorithm.  arXiv:1411.0095.
  * Chang, Sundararajan, and Keerthi (2013) Chang, K.-W.; Sundararajan, S.; and Keerthi, S. S. 2013.  Tractable semi-supervised learning of complex structured prediction models.  In _Joint European Conference on Machine Learning and Knowledge Discovery in Databases_ , 176–191. Springer.
  * Cormen et al. (2001) Cormen, T. H.; Leiserson, C. E.; Rivest, R. L.; and Stein, C. 2001.  _Introduction to Algorithms_.  The MIT Press, 2 edition.
  * Dai and Fletcher (2006) Dai, Y.; and Fletcher, R. 2006.  New algorithms for singly linearly constrained quadratic programs subject to lower and upper bounds.  _Math. Program._ , 106(3): 403–421.
  * Danskin (2012) Danskin, J. M. 2012.  _The theory of max-min and its application to weapons allocation problems_ , volume 5.  Springer Science & Business Media.
  * Duchi et al. (2008) Duchi, J.; Shalev-Shwartz, S.; Singer, Y.; and Chandra, T. 2008.  Efficient projections onto the l1-ball for learning in high dimensions.  In _Proceedings of the 25th international conference on Machine learning_ , 272–279.
  * Edelman, Ostrovsky, and Schwarz (2007) Edelman, B.; Ostrovsky, M.; and Schwarz, M. 2007.  Internet advertising and the generalized second-price auction: Selling billions of dollars worth of keywords.  _American economic review_ , 97(1): 242–259.
  * Fleder and Hosanagar (2009) Fleder, D.; and Hosanagar, K. 2009.  Blockbuster culture’s next rise or fall: The impact of recommender systems on sales diversity.  _Management science_ , 55(5): 697–712.
  * Friedlander and Tseng (2008) Friedlander, M. P.; and Tseng, P. 2008.  Exact regularization of convex programs.  _SIAM Journal on Optimization_ , 18(4): 1326–1350.
  * Gilbert (1966) Gilbert, E. 1966.  Minimizing the quadratic form on a convex set.  _SIAM J. Control_ , 4: 61–79.
  * (19) Gilbert, E.; Johnson, D.; and Keerthi, S. ????  A Fast Procedure for Computing the Distance Between Complex Objects in Three-Dimensional Space.  _IEEE Trans. Rob & Aut._, 1.
  * Grünbaum et al. (2003) Grünbaum, B.; Kaibel, V.; Klee, V.; and Ziegler, G. M. 2003.  _Convex polytopes_.  New York: Springer.  ISBN 9780387004242.
  * Guillemin and Pollack (1974) Guillemin, V.; and Pollack, A. 1974.  _Differential Topology_.  Prentice-Hall.
  * Harper and Konstan (2015) Harper, F. M.; and Konstan, J. A. 2015.  The movielens datasets: History and context.  _Acm transactions on interactive intelligent systems (tiis)_ , 5(4): 1–19.
  * Haubold et al. (2018) Haubold, C.; Uhlmann, V.; Unser, M.; and Hamprecht, F. A. 2018.  Diverse M-Best Solutions by Dynamic Programming.  _CoRR_ , abs/1803.05748.
  * He, Leng, and Li (2006) He, B.; Leng, G.; and Li, K. 2006.  Projection problems for symmetric polytopes.  _Advances in Mathematics_ , 207(1): 73–90.
  * Keerthi, Sellamanickam, and Shevade (2012) Keerthi, S. S.; Sellamanickam, S.; and Shevade, S. 2012.  Extension of TSVM to Multi-Class and Hierarchical Text Classification Problems With General Losses.  In _Proceedings of COLING 2012: Posters_ , 1091–1100. Mumbai, India: The COLING 2012 Organizing Committee.
  * Keerthi and Tomlin (2007) Keerthi, S. S.; and Tomlin, J. A. 2007.  Constructing a maximum utility slate of on-line advertisements.  _CoRR_ , abs/0706.1318.
  * Kuhn (1955) Kuhn, H. W. 1955.  The Hungarian method for the assignment problem.  _Naval research logistics quarterly_ , 2(1-2): 83–97.
  * Lin et al. (2019) Lin, Y.; Xia, Q.; He, W.; and Zhang, Q. 2019.  A Fast Iterative Check Polytope Projection Algorithm for ADMM Decoding of LDPC Codes by Bisection Method.  _IEICE Transactions on Fundamentals_ , E102–A(10).
  * Makari et al. (2013) Makari, F.; Awerbuch, B.; Gemulla, R.; Khandekar, R.; Mestre, J.; and Sozio, M. 2013\.  A distributed algorithm for large-scale generalized matching.
  * Mangasarian and Meyer (1979) Mangasarian, O. L.; and Meyer, R. 1979.  Nonlinear perturbation of linear programs.  _SIAM Journal on Control and Optimization_ , 17(6): 745–752.
  * Martins et al. (2015) Martins, A. F.; Figueiredo, M. A.; Aguiar, P. M.; Smith, N. A.; and Xing, E. P. 2015\.  Ad3: Alternating directions dual decomposition for map inference in graphical models.  _The Journal of Machine Learning Research_ , 16(1): 495–545.
  * Meyer and Polak (1970) Meyer, G.; and Polak, E. 1970.  A decomposition algorithm for solving a class of optimal control problems.  _J. Math. Anal. Apple._ , 32: 118–140.
  * Mitchell, Dem’yanov, and Malozemov (1974) Mitchell, B.; Dem’yanov, V.; and Malozemov, V. 1974.  Finding the point of a polyhedron closest to the origin.  _SIAM J. Control_ , 12: 19–26.
  * Murty (1968) Murty, K. G. 1968.  An Algorithm for Ranking All the Assignment in Order of Increasing Cost.  _Operations Research_ , 16.
  * O’Donoghue et al. (2019) O’Donoghue, B.; Chu, E.; Parikh, N.; and Boyd, S. 2019.  SCS: Splitting Conic Solver, version 2.1.2.  https://github.com/cvxgrp/scs.
  * Parikh and Boyd (2014) Parikh, N.; and Boyd, S. 2014.  Block splitting for distributed optimization.  _Mathematical Programming Computation_ , 6(1): 77–102.
  * Ramanath et al. (2021) Ramanath, R.; Keerthi, S. S.; Pan, Y.; Salomatin, K.; and Basu, K. 2021.  DuaLip: Dual Decomposition based Linear Program Solver.  https://github.com/linkedin/DuaLip/blob/master/DuaLip_AAAI_fullPaper.pdf.
  * Rambau and Ziegler (1996) Rambau, J.; and Ziegler, G. M. 1996.  Projections of polytopes and the generalized Baues conjecture.  _Discrete & Computational Geometry_, 16(3): 215–237.
  * Rush and Collins (2012) Rush, A. M.; and Collins, M. 2012.  A tutorial on dual decomposition and lagrangian relaxation for inference in natural language processing.  _Journal of Artificial Intelligence Research_ , 45: 305–362.
  * Smith (2011) Smith, N. A. 2011.  _Linguistic Structure Prediction_.  Synthesis Lectures on Human Language Technologies. Morgan and Claypool.
  * van den Berg and Friedlander (2008) van den Berg, E.; and Friedlander, M. 2008.  Probing the Pareto frontier for basis pursuit solutions.  _SIAM J. Sci. Comput._ , 31: 890––912.
  * Wang et al. (2013) Wang, J.; Zhou, J.; Wonka, P.; and Ye, J. 2013.  Lasso Screening Rules via Dual Polytope Projection.  In _NIPS_ , 1070–1078. Citeseer.
  * Wei and Banihashemi (2017) Wei, H.; and Banihashemi, A. H. 2017.  An iterative check polytope projection algorithm for ADMM-based LP decoding of LDPC codes.  _IEEE Communications Letters_ , 22(1): 29–32.
  * Wolfe (1976) Wolfe, P. 1976.  Finding the Nearest Point in A Polytope.  _Math. Program._ , 11(1): 128–149.

##  Appendix A Proofs

This section contains various details related to §3 and §4 of the paper.

A result is said to hold in general position if the result holds almost surely
in the space of perturbations and other parameters; in this section, these are
$\delta c$, $\delta b$, $\delta d$, $\delta e$, $\delta p$, $\lambda$ and $z$.
We will simply use the phrase ‘general position’ and it will be clear from the
context as to which perturbations and parameters it means.

### Proof of part 1 of Theorem 1

The proof is based on the following lemma, which we establish later.

Lemma B.1 Consider the LP

$\min_{x}~{}~{}(c+\delta c)^{T}x~{}\mbox{ s.t. }~{}Px\leq p+\delta
p,~{}Ex=e+\delta e$ (14)

where the feasible set is a compact subset of $R^{n}$ and ${\cal{P}}$ is any
open set of allowed perturbation values $(\delta c,\delta p,\delta e)$. Then,
in general position the solution $x$ of (14) is a unique point at which
exactly $n$ constraints of $Px\leq p+\delta p$, $Ex=e+\delta e$ are active. ∎

We begin with a clarification concerning the mean optimal corral dimension
associated with $x^{*}=\arg\min_{x}~{}c^{T}x~{}\mbox{s.t.}~{}Ax\leq
b,~{}x_{i}\in{\mathcal{C}}_{i},i\in[I]$. Though $x^{*}$ does not originate
from a projection operation, it is actually a projection. This comes from
Lemma 2 of [5], which implies that there exists $\bar{\gamma}>0$ such that for
all $0<\gamma\leq\bar{\gamma}$, $x^{*}$ is the projection associated with the
determination of $g_{\gamma}(\lambda_{\gamma})$ where
$\lambda_{\gamma}=\arg\max_{\lambda}g_{\gamma}(\lambda)$. Thus, it is sensible
to talk of the mean optimal corral dimension, $\mu(x^{*})$.

To prove part 1 of Theorem 1, let us apply Lemma B.1 to the case where we
combine $Ax\leq b+\delta b$ and $Dx\leq d+\delta d$ to form $Px\leq p+\delta
p$. The equality constraint, $Ex=e+\delta e$ stays as it is. The lemma implies
that in general position the solution $x^{*}$ of the LP is unique and the
number of active constraints is exactly equal to $n=IK$. Among the active
constraints, some of them, call this number $m_{act}$, come from $Ax\leq b$.
The rest of the active constraints come from $Dx\leq d$, and $Ex=e+\delta e$
i.e., the ${\mathcal{C}}_{i}$, which, in terms of optimal corral dimensions,
is equal to555The number of active constraints is related to the
representation of ${\mathcal{C}}_{i}$ by inequality constraints while optimal
corral dimension is related to the representation using vertices. We are able
to connect the two due to the special nature of the ${\mathcal{C}}_{i}$. For
the General Polytope case we would need $Dx\leq e+\delta e$ and and
$Ex=e+\delta e$ to be a minimal representation of ${\mathcal{C}}$.
$\sum_{i}(K-\mu_{i})$ for Simplex-Eq (Iq), Assignment and Slots, and upper
bounded by $\sum_{i}(K-\mu_{i})$ for Box-Cut-Eq (Iq). (It is easy to see that
these structures are unchanged if the perturbation $(\delta d,\delta e)$ is
small.) Thus,

$n\leq m_{act}+\sum_{i}(K-\mu_{i})=m_{act}+n-n\mu(x)$ (15)

Since $m_{act}\leq m$, we get $\mu(x^{*})\leq m/n$.

### Proof of part 2 of Theorem 1

Similar to Lemma B.1 we can prove the following lemma.

Lemma B.2 Consider the LP

$\min_{x}~{}~{}(c+A^{T}\lambda+\delta c)^{T}x~{}\mbox{ s.t. }~{}Dx\leq
d+\delta d,~{}Ex=e+\delta e$ (16)

where the feasible set is a compact subset of $R^{n}$ and ${\cal{D}}$ is any
open set of allowed perturbation/parameter values $(\delta c,\delta d,\delta
e,\lambda)$. Then, in general position the solution $x$ of (16) is a unique
point at which exactly $n$ constraints of $Dx\leq d+\delta d$, $Ex=e+\delta e$
are active. ∎

To prove part 2 of Theorem 1 we just apply Lemma B.2. Since $Ax\leq b$
constraint is absent, analysis similar to the proof of part 1 yields that in
general position we have $\mu(x_{0}^{*}(\lambda))=0$.

### Proofs of Lemmas B.1 and B.2

The two lemmas are geometrically intuitive; however, proving them rigorously
is somewhat technical and detailed.

Let us start by proving Lemma B.1. The Lagrangian for the LP in (14) with dual
vectors $\psi$ and $\phi$ is

$L=(c+\delta c)^{T}x+\psi^{T}(Px-p-\delta p)+\phi^{T}(Ex-e-\delta e)$ (17)

A triple $(x,\psi,\phi)$ is optimal iff the following KKT conditions hold:

$\displaystyle P^{T}\psi+E^{T}\phi+c+\delta c=0,~{}Px\leq p+\delta p,$ (18)
$\displaystyle\psi\geq 0,~{}\psi^{T}(Px-p-\delta p)=0,~{}Ex=e+\delta e$ (19)

Let

$S_{1}=\\{r:(Px-p-\delta p)_{r}=0\\},~{}S_{2}=\\{r\in S_{1}:\psi_{r}=0\\}$
(20)

$S_{1}$ is the set of active (primal) inequality constraints and $S_{2}\subset
S_{1}$ is the subset with positive $\psi$ values. Note that not all the
$\psi_{r},r\in S_{1}$ are required to be positive. Note also that
$\psi_{r}=0\;\forall r\not\in S_{1}$. Let us rewrite the KKT conditions in
(18) and (19) as

$\displaystyle P^{T}\psi+E^{T}\phi+c+\delta c=0,~{}(Px-p-\delta p)_{S_{1}}=0,$
(21) $\displaystyle(Px-p-\delta p)_{S_{1}^{c}}<0,~{}(\psi)_{S_{1}\backslash
S_{2}}>0,~{}(\psi)_{S_{1}^{c}\cup S_{2}}=0$ (22) $\displaystyle Ex-e-\delta
e=0$ (23)

In the above, we are analyzing what happens at a solution $(x,\psi,\phi)$. Let
us turn this around a bit and look at the (finite) collection of all
$(S_{1},S_{2})$ possibilities satisfying $S_{2}\subset S_{1}$, and, for a
given $(S_{1},S_{2})$ look at the set of all optimal solutions $(x,\psi,\phi)$
satisfying (21)-(23), call it ${\cal{O}}(S_{1},S_{2})$. Note that
${\cal{O}}(S_{1},S_{2})$ will turn out to be empty for most $(S_{1},S_{2})$
pairs, which is fine. Using this definition, the set of all possible solutions
of the LP in (14) can be written as

${\cal{O}}=\cup_{(S_{1},S_{2})}{\cal{O}}(S_{1},S_{2})$ (24)

Let us take one $(S_{1},S_{2})$. We will use Transversality theorem (a
corollary of Sard’s theorem) to show that ${\cal{O}}(S_{1},S_{2})$ is at most
a singleton in general position. We state Transversality theorem [21] in the
simplest form that suffices our purpose.

Transversality theorem Suppose $F:Y\times Z\to R^{N}$ is a smooth
differentiable function such that (a) $Y$ and $Z$ are open sets in $R^{L}$ and
$R^{M}$ and (b) the Jacobian of $F$ with respect to $(y,z)$ has full row rank
for all $(y,z)\in Y\times Z$. Then in general position the set
$\tilde{Y}(z)=\\{y\in Y:F(y,z)=0\\}$ is either empty or a differentiable
manifold of dimension $(L-N)$. ∎

To set up Transversality theorem for proving Lemma B.1, let $y=(x,\psi,\phi)$,
$z=(\delta c,\delta p,\delta e)$, $Z={\cal{P}}$,

$\displaystyle Y=\\{(x,\psi,\phi):(Px-p-\delta
p)_{S_{1}^{c}}<0,~{}(\psi)_{S_{1}\backslash S_{2}}>0\\}$ (25)

and the function, $F:Y\times Z\to R^{N}$ as

$F=\left(\begin{array}[]{c}P^{T}\psi+E^{T}\phi+c+\delta c\\\ (Px-p-\delta
p)_{S_{1}}\\\ Ex-e-\delta e\\\ (\psi)_{S_{1}^{c}\cup S_{2}}\end{array}\right)$
(26)

so that ${\cal{O}}(S_{1},S_{2})$ can be written as

${\cal{O}}(S_{1},S_{2})=\tilde{Y}(z)=\\{y\in Y:F(y,z)=0\\}$ (27)

Note that $Y$ and $Z$ are open sets. If we look at the Jacobian of $F$ with
respect to the subset of variables, $(\delta c,(\delta p)_{S_{1}},\delta
e,(\psi)_{S_{1}^{c}\cup S_{2}})$ it is a non-singular diagonal matrix with
$1$’s and $-1$’s in the diagonal. Hence the Jacobian of $F$ with respect to
$(y,z)$ has full row rank. Thus, by Transversality theorem, in general
position the set ${\cal{O}}(S_{1},S_{2})$ is either empty or a differentiable
manifold of dimension $(L-N)$. With $n$ being the number of variables in $x$
and $\ell$ being the number of constraints in $Px\leq p+\delta p$, $Ex\leq
e+\delta e$ note that $L=n+\ell,~{}N=n+\ell+|S_{2}|$. Let’s consider two
cases.

Case 1. $\mathbf{|S_{2}|>0.}$ Since $(L-N)$ is negative, it automatically
means that in general position the set ${\cal{O}}(S_{1},S_{2})$ is empty.
Thus, in general position, we don’t have to consider the case
$S_{2}\not=\emptyset$.

Case 2. $\mathbf{|S_{2}|=0.}$ For this case, $(L-N)=0$ and, in general
position the set ${\cal{O}}(S_{1},S_{2})$ is a zero dimensional manifold,
i.e., a set of isolated points. By convexity of ${\cal{O}}(S_{1},S_{2})$, it
can be at most a singleton. Also, if ${\cal{O}}(S_{1},S_{2})\not=\emptyset$,
then, since $(Px-p-\delta p)_{S_{1}}=0$ and $Ex-e-\delta e=0$ independently
determine $x$, and, $P^{T}\psi+c+\delta c=0$, $(\psi)_{S_{1}^{c}}=0$ together
determine $\psi$, we must have $|S_{1}|=n$; otherwise, in general position,
one of the two equation sets will have no solution.

Let us now come to ${\cal{O}}$ in (24). In general position, convexity of
${\cal{O}}$ implies that the only possibility is that there is a unique
$S_{1}$ such that
${\cal{O}}=O(S_{1},\emptyset)=\\{(x^{*},\psi^{*},\phi^{*})\\}$, the unique
optimal solution of (14). ∎

The lines of proof for Lemma B.2 are similar. $Px\leq p+\delta p$ is replaced
by $Dx\leq e+\delta e$ and there is an additional parameter $\lambda$ that is
included in $z$.

### Wolfe’s algorithm

In §3 we described the essence of Wolfe’s algorithm. Here we give fuller
details in Algorithm 3.

Steps 5-22 form one iteration of the algorithm going from one corral to a
better corral. Steps 5 and 6 correspond to checking optimality of $x$, the
projection to the current best corral; these steps are same as each of (13)
and (12). Step 5 forms the most expensive part of the algorithm. Note that, in
problems with structure in the vertices of ${\mathcal{C}}$, step 5 can be done
very efficiently without enumerating all vertices. As we explained in §4, this
is exactly what happens in Assignment , Slate and Box-Cut-Eq (Iq).

If checking of optimality of $x$ fails, the algorithm enters the minor cycle
in which a next, better corral is identified efficiently. Steps 11-22 form the
minor cycle. Step 12 requires the projection of $\hat{x}$ to $\tt{aff}(S)$,
the affine hull of the vertices in $S$, which can be obtained by solving a
system of linear equations in $(|S|+1)$ unknowns. Since the corral dimensions
$(|S|-1)$ encountered in marketplace LPs are small, the affine hull
computation is cheap.

In many cases of ${\mathcal{C}}$, the minor cycle breaks at the first
cycle/entrance to step 14. The modified Duchi et al algorithm for Simplex-Eq
(algorithm 2) is an example. This is also true in most iterations of other
cases of ${\mathcal{C}}$.

Algorithm 3 Wolfe’s algorithm

1:Input:  $\hat{x}$, $\mathcal{B}=$ Vertex set of ${\mathcal{C}}$

2:Output:  $x^{*}=\arg\min_{x\in{\mathcal{C}}}\|x-\hat{x}\|^{2}$

3:In all iterations we maintain $S=\\{v_{r}\\}$, a subset of vertices, and
$x=\sum_{r\in S}\rho_{r}v_{r}$ as a convex combination of the vertices in $S$.
At the end $S$ forms the optimal corral.

4:Find $v$ the vertex closest to $\hat{x}$.

5:Initialize $x\leftarrow v$, $S\leftarrow\\{v\\}$ and $\rho_{1}\leftarrow 1$.

6:while true do$\triangleright$ Major Cycle

7: $v\mathrel{:}=\arg\min_{p\in\mathcal{B}}(x-\hat{x})^{T}p$

8: if $(x-\hat{x})^{T}x\leq(x-\hat{x})^{T}v$ then

9: Set $x^{*}=x$.

10: break

11: end if

12: $S\mathrel{:}=S\cup\\{v\\}$, $\rho_{|S|}\leftarrow 0$

13: while true do$\triangleright$ Minor Cycle

14: Find $y=\arg\min_{z\in\tt{aff}(S)}\|z-\hat{x}\|^{2}$ and

15: $\alpha$ s.t. $y=\sum_{r\in S}\alpha_{r}v_{r}$

16: if $\alpha_{r}\geq 0\;\forall_{r}$ then

17: break$\triangleright$ $y\in\tt{co}(S)$, so end minor cycle

18: else$\triangleright$ If $y\not\in\tt{conv}(S)$, then update $x$ to the
intersection of the boundary of $\tt{co}(S)$ and the segment joining $y$ and
previous $x$. Delete points from $S$ which are not required to describe the
new $x$ as a convex combination.

19:
$\theta\mathrel{:}=\min_{r:\alpha<0}\frac{\rho_{r}}{\rho_{r}-\alpha_{r}}$$\triangleright$
Since $x=\sum_{r}\rho_{r}v_{r}$

20: $x\leftarrow\theta y+(1-\theta)x$$\triangleright$ Using $\theta$, the new
$x$ lies in $\tt{conv}(S)$.

21: $\rho_{r}\leftarrow\theta\alpha_{r}+(1-\theta)\rho_{r}\;\forall
r$$\triangleright$ Sets the coefficients of the new $x$

22: $S\leftarrow\\{r:\rho_{r}>0\\}$$\triangleright$ Delete points which have
$\rho_{r}=0$. This deletes at least one point.

23: end if

24: end while$\triangleright$ After the minor loop terminates, $x$ is updated
to be the affine minimizer of the current set $S$ $x\leftarrow y$

25:end while

26:return $x^{*}$