Faster width-dependent algorithm for mixed packing and covering LPs
Digvijay Boob Saurabh Sawlani Di Wang






Faster Width-dependent Algorithm for Mixed Packing and Covering LPs
    
	Digvijay Boob 

	Georgia Tech

	Atlanta, GA 

	

	Saurabh Sawlani 

	Georgia Tech

	Atlanta, GA 

	
	Di WangWork done when author was at Georgia Tech. 

	Google AI

	Atlanta, GA 

	

    December 30, 2023
============================================================================================================================================================================================






	In this paper, we give a faster width-dependent algorithm for mixed packing-covering LPs.
	Mixed packing-covering LPs are fundamental to combinatorial optimization in computer science and operations research.
	
	Our algorithm finds a 1+ approximate solution in time O(Nw/ ), where N is number of nonzero entries in the constraint matrix,
	and w is the maximum number of nonzeros in any constraint.
	This run-time is better than Nesterov's smoothing algorithm which requires O(N√(n)w/ ) where n is the dimension of the problem.
	Our work utilizes the framework of area convexity introduced in [Sherman-FOCS'17]
	to obtain the best dependence on 
	while breaking the infamous ℓ_∞ barrier to eliminate the factor of √(n).
	The current best width-independent algorithm for this problem runs in time O(N/^2) [Young-arXiv-14]
	and hence has worse running time dependence on .
	Many real life instances of the mixed packing-covering problems exhibit small width and for such cases, our algorithm can report higher precision results when compared to width-independent algorithms.
	As a special case of our result, we report a 1+ approximation algorithm for the densest subgraph problem which runs in time O(md/ ), where m is the number of edges in the graph and d is the maximum graph degree.








arabic





§ INTRODUCTION
 

Mixed packing and covering linear programs (LPs) are
a natural class of LPs
where coefficients, variables, and constraints are
non-negative.
They model a wide range of important problems
in combinatorial optimization and
operations research.
In general, they model any problem
which contains a limited set of available
resources (packing constraints)
and a set of demands to fulfill (covering constraints).

Two special cases of the problem
have been widely studied in literature: pure packing,
formulated as max_x {b^T x | Px ≤ p };
and pure covering, formulated as
min_x {b^T x | Cx ≥ c } where P,p,C, c, b
are all non-negative.
These are known to model fundamental problems
such as maximum bipartite graph matching, minimum set cover,  etc <cit.>.
Algorithms to solve packing and covering LPs
have also been applied to great effect in designing
flow control systems <cit.>,
scheduling problems <cit.>, zero-sum matrix games <cit.>
and in mechanism design <cit.>.
In this paper, we study the mixed packing and covering (MPC) problem, formulated as checking the feasibility of the set:
{x | Px ≤ p, Cx ≥ c},
where P,C,p,c are non-negative.
We say that x is an -approximate solution to MPC
if it belongs to the relaxed set
{x | Px ≤ (1+)p, Cx ≥ (1-)c}.
MPC is a generalization of pure packing and pure covering,
hence it is applicable to a wider range of problems
such as multi-commodity flow on graphs <cit.>,
non-negative linear systems and X-ray tomography <cit.>.



General LP solving techniques such as the interior point method
can approximate solutions to MPC
in as few as O(log(1/)) iterations - however,
they incur a large per-iteration cost.
In contrast, iterative approximation algorithms
based on first-order optimization methods require poly(1/) iterations,
but the iterations are fast and in most cases are conducive to efficient parallelization.
This property is of utmost importance in the context
of ever-growing datasets and the availability of
powerful parallel computers, resulting in much faster algorithms
in relatively low-precision regimes.



 §.§ Previous work


In literature, algorithms for the MPC problem
can be grouped into two broad categories:
width-dependent and width-independent.
Here, width is an intrinsic property of a linear program
which typically
depends on the dimensions and the largest entry of the constraint matrix,
and is an indication of the range of values any constraint
can take.

In the context of this paper and the MPC problem, we define w_P and w_C as the maximum
number of non-zeros in any constraint in P and C respectively.
We define the width of the LP as w_maxmax(w_P, w_C).

One of the first approaches used to solve LPs was Langrangian-relaxation:
replacing hard constraints with loss functions
which enforce the same constraints indirectly.
Using this approach, Plotkin, Schmoys and Tardos <cit.>
and Grigoriadis and Khachiyan <cit.>
obtained width-dependent polynomial-time approximation algorithms for MPC.
Luby and Nisan <cit.> gave the first
width-dependent parallelizable algorithm for pure packing and pure covering,
which ran in O(^-4) parallel time,
and O(N^-4) total work.
Here, parallel time (sometimes termed as depth) refers to the longest
chain of dependent operations, and work refers to the
total number of operations in the algorithm.


Young <cit.> extended this technique to give
the first width-independent parallel algorithm
for MPC in O(^-4) parallel time,
and O(m d ^-2) total work[d here
is the maximum number of constraints that any variable appears in.].
Young <cit.> later improved his algorithm
to run using total work O(N ^-2).
Mahoney <cit.> later gave an algorithm with a faster parallel
run-time of O(^-3).

The other most prominent approach in literature towards solving an LP is by converting it into a smooth function <cit.>,
and then applying general first-order optimization techniques <cit.>.
Although the dependence on  from using first-order
techniques is much improved,
it usually comes at the cost of sub-optimal dependence on the input size and width.
For the MPC problem,
Nesterov's accelerated method <cit.>, as well as Bienstock and Iyengar's adaptation <cit.> of Nesterov's smoothing <cit.>,
give rise to algorithms with runtime linearly depending on ^-1,
but with far from optimal dependence on input size and width.
For pure packing and pure covering problems, however,
Allen-Zhu and Orrechia <cit.> were the first to incorporate Nesterov-like acceleration while still being able to obtain near-linear width-independent runtimes, giving a O(N^-1) time algorithm for the packing problem.
For the covering problem, they gave a
O(N^-1.5) time algorithm, which was then improved to O(N^-1) by <cit.>.
Importantly, however, the above algorithms do not generalize to MPC.














 §.§ Our contributions


We give the best parallel width-dependent algorithm for MPC,
while only incurring a linear dependence on ^-1
in the parallel runtime and total work.
Additionally, the total work has near-linear dependence on
the input-size.
Formally, we state our main theorem as follows.

 
There exists a parallel -approximation algorithm for the mixed packing covering
problem, which runs in O(w ·^-1) parallel time,
while performing O(w · N ·^-1) total work,
where N is the total number of non-zeros in the constraint matrices,
and w is the width of the given LP.


Table <ref> compares the running time of our algorithm
to previous works solving this problem (or its special cases).





Sacrificing width independence for faster convergence with respect to precision
proves to be a valuable trade-off for several combinatorial optimization
problems which naturally have a low width.
Prominent examples of such problems which are not pure packing or covering
problems include multicommodity flow and densest subgraph,
where the width is bounded by the degree of a vertex.
In a large number of real-world graphs, the maximum vertex degree is usually
small, hence our algorithm proves to be much faster when we want high-precision
solutions. We explicitly show that this result directly gives
the fastest algorithm for the densest subgraph problem on low-degree graphs in Appendix <ref>.


§ NOTATION AND DEFINITIONS
 
For any integer q, we represent using ·q
the q-norm of any vector.
We represent the infinity-norm as ·∞.
We denote the infinity-norm ball (sometimes called the ℓ_∞ ball) as the set _∞^n(r) { x ∈^n: x∞≤ r }. The nonnegative part of this ball is denoted as _+,∞^n(r) = {x ∈^n: x ≥_n, x∞≤ r }. For radius r = 1, we drop the radius specification and use a short notation _∞^n and _+,∞^n.  We denote the extended simplex of dimension k as Δ_k^+ { x ∈^k: ∑_i = 1^k x_i ≤ 1 }. For any y ≥_k, _Δ_k^+(y) = y / y1 if y1≥ 1.

Further, for any set K, we represent its interior, relative interior and closure as (K), (K) and (K), respectively. Function exp is applied to a vector element wise. Division of two vectors of same dimension is also performed element wise.

For any matrix A, we use (A) to denote the number of nonzero entries in it.
We use A_i,: and A_:, j to refer to the ith row and jth column of A respectively. We use notation A_ji or A_j,i alternatively to denote element in j-th row and i-th column of matrix A.
A∞ denotes the operator norm A∞→∞sup_x ≠ 0Ax∞/x∞.
For a symmetric matrix A and an antisymmetric matrix B, we define an operator  as 
A  B ⇔[  A -B;  B  A ] is positive semi-definite.

We formally define an -approximate solution to the mixed packing-covering
(MPC) problem as follows.

	We say that x is an -approximate solution of the mixed packing-covering problem if x satisfies x ∈_+,∞^n(1), Px ≤ (1+) _p and Cx ≥ (1-)_c.

Here, _k denotes a vectors of 1's of dimension k for any integer k.

The saddle point problem on two sets x ∈ X and y ∈ Y can be defined as follows:

    min_x ∈ Xmax_y ∈ Y (x,y)

where (x,y) is some bilinear form between x and y.
For this problem, we define the primal-dual gap function
as sup_(x, y) ∈ X × Y(x, y)  - (x, y).
This gap function can be used as measure of accuracy of the above saddle point solution.

We say that  (x,y) ∈ X × Y is an -optimal solution for (<ref>) if sup_(x, y) ∈ X × Y(x, y)  - (x, y) ≤.




§ TECHNICAL OVERVIEW
 
The mixed packing-covering (MPC) problem is formally defined as follows.

Given two nonnegative matrices P ∈^p × n, C ∈^c × n, find an x ∈^n, x ≥, x∞≤ 1 such that Px ≤_p and Cx ≥_c if it exists, otherwise report infeasibility.

Note that the vector of 1's on the right hand side of packing and covering constraints
can be obtained by simply scaling each constraint appropriately.
We also assume that each entry in the matrices P and C is at most one.
This assumption, and subsequently the ℓ_∞ constraints on x also cause no loss of generality.[This transformation can be achieved by adapting techniques from <cit.> while increasing dimension of the problem up to a logarithmic factor.
Details of this fact are in the Appendix <ref>
in the full paper (supplementary file).
For the purpose of the main text, we work with this assumption.]

We reformulate MPC as a saddle point problem, as defined in Section <ref>;

    λ^* min_x ∈_+,∞^n  max_y ∈Δ^+_c, z ∈Δ^+_p  L(x,y,z),

where  L(x,y,z) := [ y^T  z^T] [   P -_p;  -C  _c ][ x; 1 ].
 The relation between the two formulation is shown in Section <ref>.
 



For the rest of the paper, we focus on the saddle point formulation (<ref>).
 η(x) max_y ∈Δ^+_c,  z ∈Δ^+_p L(x,y,z) is a piecewise linear convex function.
Assuming oracle access to this “inner" maximization problem,
the “outer" problem of minimizing η(x) can be performed using
first order methods like mirror descent, which are suitable
when the underlying problem space is the unit ℓ_∞ ball.
One drawback of this class of methods is that their rate of convergence,
which is standard for non-accelerated first order methods on non-differentiable objectives,
is O(1/^2) to obtain an -approximate minimizer x
of η which satisfies η(x) ≤η^* +, where η^* is the optimal value.
This means that the algorithm needs to access the inner maximization oracle
O(1/^2) times,
which can become prohibitively large in the high precision regime.

Note that even though η is a piecewise linear non-differentiable function, it is not a black box function, but a maximization linear functions in x.
This structure can be exploited using Nesterov's smoothing technique <cit.>.
In particular, η(x) can be approximated by choosing a strongly convex[3] function ϕ: Δ_p^+ ×Δ_c^+ → and considering

    η (x)= max_y ∈Δ^+_c,  z ∈Δ^+_p  L(x,y,z) - ϕ(y,z).

This strongly convex regularization yields that η is a Lipschitz-smooth[Definitions of Lipschitz-smoothness and strong convexity can be found in many texts in nonlinear programming and machine learning. e.g. <cit.>. Intuitively, f is Lipschitz-smooth if the rate of change of ∇ f can be bounded by a quantity known as the “constant of Lipschitz smoothness”.] convex function.
If L is the constant of Lipschitz smoothness of η then application of any of the accelerated gradient methods in literature will converge in O(√(L/)) iterations.
Moreover, it can also be shown that in order to construct a smooth -approximation η of η,
the Lipschitz smoothness constant L can be chosen to be of the order O(1/),
which in turn implies an overall convergence rate of O(1/).
In particular, Nesterov's smoothing achieves an oracle complexity of O((P∞+C∞)  D_x max{D_y, D_z}^-1 ), where
where D_x, D_y and D_z denote the sizes of the ranges of their respective
regularizers which are strongly convex functions.
D_y and D_z can be made of the order of logp and logc, respectively.
However, D_x can be problematic since x belongs to an ℓ_∞ ball.
More on this will soon follow.

Nesterov's dual extrapolation algorithm<cit.> gives a very similar complexity but
is a different algorithm in that it directly addresses the saddle point formulation (<ref>) rather than viewing the problem as
optimizing a non-smooth function η.
The final convergence for the dual extrapolation algorithm is given in terms of the primal-dual gap function of the saddle point problem (<ref>).
This algorithms views the saddle point problem as solving variational inequality for an appropriate monotone operator in joint domain (x,y,z). 
Moreover, as opposed to smoothing techniques which only regularize the dual, this algorithm regularizes both primal and dual parts, hence is a different scheme altogether.

Note that for both schemes mentioned above, the maximization oracle itself has an analytical expression which involves matrix-vector multiplication.
Hence each call to the oracle incurs a sequential run-time of (P)+(C).
Then, overall complexity for both schemes is of order O(((P)+(C))(P∞+C∞)  D_xmax{D_y, D_z}^-1).


  §.§.§ The ℓ_∞ barrier

Note that the  both methods, i.e.,
Nesterov's smoothing and dual extrapolation, involves a D_x term, which denotes the range of a convex function over the domain of x.
The following lemma states a lower bound for this range in case of ℓ_∞ balls.


Any strongly convex function has a range of at least Ω(√(n))
on any ℓ_∞ ball.





Since D_x ≥√(n)/2√(2) for each member function of this wide class,
there is no hope of eliminating this √(n) factor using techniques
involving explicit use of strong convexity.

So, the goal now is to find a function with a small range over ℓ_∞ balls,
but still act as good enough regularizers to enable accelerated convergence of
the descent algorithm.
In pursuit of breaking this ℓ_∞ barrier, we draw inspiration
from the notion of area convexity introduced by Sherman <cit.>.
Area convexity is a weaker notion than strong convexity, however, it
is still strong enough to ensure that accelerated first order methods still go through when
using area convex regularizers. Since this is a weaker notion than strong convexity, we can construct area convex functions which have range of O(n^o(1)) on ℓ_∞ ball.

First, we define area convexity, and then go on to mention its relevance
to the saddle point problem (<ref>).

Area convexity is a notion defined in context of a matrix A ∈^a × b and a convex set K ⊆^a+b. Let M_A [ _b × b   -A^T;      A _a × a ].


 	A function ϕ is area convex with respect to a matrix A on a convex set K iff for any t,u,v ∈ K, ϕ satisfies ϕ[]t+u+v3≤13( ϕ(t) +ϕ(u)+ϕ(v) ) -13√(3)(v-u)^TM_A(u-t)
 
 To understand the definition above, first lets look at the notion of strong convexity.
 ϕ is strongly convex if for any two points t,u, 1/2(ϕ(t)+ ϕ(u)) exceeds ϕ1/2(t+u) by an amount proportional to t-u2^2.
 Definition <ref> generalizes this notion in context of matrix A for any three points x, y, z.
 ϕ is area-convex on set K if for any three points t, u, v ∈ K, we have 1/3(ϕ(t) + ϕ(u)+ ϕ(v)) exceeds ϕ(1/3(t+u+v)) by an amount proportional to the area of the triangle defined by the convex hull of t,u,v.
 
 Consider the case that points t,u,v are collinear.
 For this case, the area term (i.e., the term involving M_A) in Definition <ref> is 0 since matrix M_A is antisymmetric.
 In this sense, area convexity is even weaker than strict convexity.
 Moreover, the notion of area is parameterized by matrix A.

 To see a specific example of this notion of area, consider
 A = [  0 -1;  1  0 ]
 and t, u,v ∈^2.
 Then, for all possible permutations of t,u,v, the area term takes a value equal to ±(t_1(u_2-v_2)+ u_1(v_2-t_2) + v_1(t_2-u_2)).
 Since the condition holds irrespective of the permutation so we must have that ϕt+u+v3≤13( ϕ(t) +ϕ(u)+ϕ(v) ) -13√(3)t_1(u_2-v_2)+ u_1(v_2-t_2) + v_1(t_2-u_2).
 But note that area of triangle formed by points t,u,v is equal to 1/2t_1(u_2-v_2)+ u_1(v_2-t_2) + v_1(t_2-u_2).
 Hence the area term is just a high dimensional matrix based generalization of the area of a triangle.
 
 Coming back to the saddle point problem (<ref>),
 we need to pick a suitable area convex function ϕ on the set
 _+,∞^n ×Δ_p^+×Δ_c^+.
 Since ϕ is defined on the joint space,
 it has the property of joint regularization vis a vis (<ref>).
 However, we need an additional parameter: a suitable matrix M_A. 
 The choice of this matrix is related to the bilinear form of the
 primal-dual gap function of (<ref>).
 We delve into the technical details of this in Section <ref>,
 however, we state that the matrix is composed of P, C and some additional constants. 
 The algorithm we state exactly follows Nesterov's dual extrapolation
 method described earlier.
 One notable difference is that in <cit.>,
 they consider joint regularization by a strongly convex function
 which does not depend on the problem matrices P,C but only on the constraint set _+,∞^n ×Δ_p^+×Δ_c^+.
 Our area convex regularizer, on the other hand, is tailor made
 for the particular problem matrices P,C as well as the constraint set.




§ AREA CONVEXITY FOR MIXED PACKING COVERING LPS
 
In this section, we present our technical results and algorithm
for the MPC problem, with the end goal of proving Theorem <ref>.
First, we relate an (1+)-approximate solution to the saddle point problem to an -approximate solution to MPC.
Next, we present some theoretical background towards the goal of choosing
and analyzing an appropriate area-convex regularizer in the context of the saddle point formulation, where the key requirement of the
area convex function is to obtain a provable and efficient convergence result.
Finally, we explicitly show an area convex function which is generated using a simple “gadget" function.
We show that this area convex function satisfies all key requirements and hence achieves the desired
accelerated rate of convergence.
This section closely follows <cit.>, in which the author chooses an area convex
function specific to the undirected multicommodity flow problem.
Due to space constraints, we relegate almost all proofs to Appendix <ref>
and simply include pointers to proofs in <cit.> when it is directly applicable. 



 §.§ Saddle Point Formulation for MPC

Consider the  saddle point formulation in (<ref>) for MPC problem.
which is closely related to mixed packing-covering problem (P). Rewrite (P) as 
    Px    ≤_p
    
    -Cx    ≤ -_c
    
    x    ∈_+,∞^n(1)

Then for any dual variables y ∈^p_+ and z ∈^c_+, we define L(x,y,z) := [ y^T  z^T] [   P -_p;  -C  _c ][ x; 1 ]. 
Then consider the following saddle point problem

    min_x ∈_+,∞^n max_y ∈Δ^+_c
      z ∈Δ^+_p  L(x,y,z)

Here _+,∞^n :={x ∈^n: x ≥, x∞≤ 1 }, Δ_p^+ := {x ∈^p: x ≥, ∑_i =1^p x_i ≤ 1} and Δ_c^+ := {x ∈^c: x ≥, ∑_i =1^c x_i ≤ 1}.
Given a feasible primal-dual feasible solution pair (x,y,z) and (x, y, z) for (<ref>),
we denote w = (x,u, y,z) and w = (x, u, y, z) where u , u∈.
Then, we define a function Q: ^n+1+p+c×^n+1+p+c→ as 

    Q(w, w)  [y^T z^T]
    [   P -_p;  -C  _c ][ x; u ]
    -
    [y^T  z^T]
    [   P -_p; - C  _c ][ x; u ].

Note that if u = u = 1, then 
    sup_w∈Q(w, w) = sup_x∈_+,∞^n, y∈Δ_p^+, z∈Δ_c^+L(x,y, z)- L(x, y, z)

is precisely the primal-dual gap function defined in Section <ref>.
Notice that if (x^*, y^*,z^*) is a saddle point of (<ref>), then we have 

    L(x^*, y , z) ≤ L(x^*,y^*, z^*) ≤ L(x, y^*,z^*)

for all x ∈_+,∞^n, y ∈Δ_p^+, z∈Δ_c^+.
From above equation, it is clear that Q(w, w^*) ≥ 0 for all w ∈ where
_+,∞^n(1)×{1}×Δ_p^+ ×Δ_c^+
and w^* = (x^*, 1, y^*, z^*) ∈.
Moreover, Q(w^*,w^*) = 0. This motivates the following accuracy measure of the candidate approximate solution w.

	We say that w ∈ is an -optimal solution of (<ref>) iff 
	
    sup_w∈Q(w, w)  ≤.


 
	Recall the definition of M_A for a matrix A in Section <ref>.
	We can rewrite Q(w,w) = w^T J w where J = M_H and
	
    H = [   P -_p;  -C  _c ]  ⇒   J := [ _n × n _n × 1   -P^T    C^T;  _1× n      0   _p^T  -_c^T;      P    -_p _p × p _p × c;     -C     _c _c × p _c × c;        ].

	Thus, the gap function in Definition <ref> can be written in the bilinear form sup_w∈w^TJw.
 


Lemma <ref> relates the -optimal solution of (<ref>) to the -approximate solution to MPC.


	Let (x,y,z) satisfy sup_(x, y, z) ∈_+,∞^n ×Δ_p^+ ×Δ_c^+ L(x, y, z) - L(x, y, z) ≤. Then either 

	1. x is an -approximate solution of MPC, or

	2. y,z satisfy y^T(Px-_p) + z^T(-Cx + _c) > 0 for all x∈_+,∞^n.

This lemma states that in order to find an -approximate solution of MPC, it suffices to find -optimal solution of (<ref>).
Henceforth, we will focus on -optimality of the saddle point formulation (<ref>).



 §.§ Area Convexity with Saddle Point Framework

Here we state some useful lemmas which help in determining whether a differentiable function is area convex.
We start with the following remark which follows from the definition of area convexity (Definition <ref>).

	If ϕ is area convex with respect to A on a convex set K,
	and K⊆ K is a convex set,
	then ϕ is area convex with respect to A on K.

The following two lemmas from <cit.> provide the key characterization of area convexity. 






	Let A ∈^2 × 2 symmetric matrix. A [  0 -1;  1  0 ]⇔ A ≽ 0 and (A) ≥ 1.

 
	Let ϕ be twice differentiable on the interior of convex set K, i.e., (K).
	
		
  * If ϕ is area convex with respect to A on (K), then d^2ϕ(x)  M_A for all x ∈(K).
		
  * If d^2ϕ(x)  M_A for all x ∈(K), then ϕ is area convex with respect to 1/3A on (K). Moreover, if ϕ is continuous on (K), then ϕ is area convex with respect to 1/3A on (K).
	

In order to handle the operator  (recall from Section <ref>),
we state some basic but important properties of this operator,
which will come in handy in later proofs.
 
	For symmetric matrices A and C and antisymmetric matrices B and D,
	
		
  * If A  B then A  (-B).
		
  * If A  B and λ≥ 0 then λ A λ B.
		
  * If A  B and C  D then A+ C   (B+D).
	


Having laid a basic foundation for area convexity,
we now focus on its relevance to solving the saddle point problem (<ref>).
Considering Remark <ref>,
we can write the gap function criterion of optimality in terms of bilinear form of the matrix J.

Suppose we have a function ϕ which is area convex with respect to H on set . 
Then, consider the following jointly-regularized version of the bilinear form:

    η(w) := sup_w∈ w^TJw - ϕ(w).

Similar to Nesterov's dual extrapolation, one can attain O(1/) convergence
of accelerated gradient descent for function η(w) in (<ref>)
over variable w.
In order to obtain gradients of η(w), we need access to
_w∈w^TJw -ϕ(w).
However, it may not be possible to find an exact maximizer in all cases.
Again, one can get around this difficulty by instead using an approximate optimization oracle
of the problem in (<ref>).

	A δ-optimal solution oracle (OSO) for ϕ:→ takes input a and outputs w ∈ such that 
	
    a^Tw -ϕ(w) ≥sup_w∈ a^Tw -ϕ(w) -δ.
 
 

Given Φ as a δ-OSO for a function ϕ,
consider the following algorithm (Algorithm <ref>):

For Algorithm <ref>, <cit.> shows the following:

	Let ϕ: →[-ρ, 0]. Suppose ϕ is area convex with respect to 2√(3)H on . Then for J = M_H and for all t≥ 1 we have w_t/ t ∈ and, 
    sup_w∈w Jw_tt≤δ + ρt.

 	In particular, in ρ/ iterations, Algorithm <ref> obtain (δ+)-solution of the saddle point problem (<ref>).

The analysis of this lemma closely follows the analysis of Nesterov's dual extrapolation.

Note that, each iteration consists of O(1) matrix-vector multiplications,
O(1) vector additions,
and O(1) calls to the approximate oracle.
Since the former two are parallelizable to O(log n) depth,
the same remains to be shown for the oracle computation
to complete the proof of the run-time in Theorem <ref>.

Recall from the discussion in Section <ref> that the
critical bottleneck of Nesterov's method is that diameter of the ℓ_∞ ball is Ω(√(n)),
which is achieved even in the Euclidean ℓ_2 norm.
This makes ρ in Lemma <ref> to also be Ω(√(n)),
which can be a major bottleneck for high dimensional LPs, which are commonplace among real-world applications.

Although, on the face of it, area convexity applied to the saddle point formulation (<ref>)
has a similar framework to Nesterov's dual extrapolation,
the challenge is to construct a ϕ for which we can overcome the above bottleneck.
Particularly, there are three key challenges to tackle:

	1. We need to show that existence of a function ϕ that is area convex with respect to H on .

	2. ϕ: → [-ρ, 0] should be such that ρ is not too large.
 
	3. There should exist an efficient δ-OSO for ϕ.

In the next subsection, we focus on these three aspects in order to complete our analysis.



 §.§ Choosing an area convex function

First, we consider a simple 2-D gadget function and prove a “nice" property of this gadget.
Using this gadget, we construct a function which can be shown to be area convex
using the aforementioned property of the gadget.

Let γ_β:_+^2 → be a function parameterized by β defined as 

    γ_β(a,b) = baloga + β blogb .





	Suppose β≥ 2. Then d^2γ_β(a,b) ≽[  0 -1;  1  0 ] for all a ∈ (0,1] and b>0.


Now, using the function γ_β,
we construct a function ϕ
and use the sufficiency criterion provided in Lemma <ref> to show that
ϕ is area convex with respect to J on .
Note that our set of interest  is not full-dimensional,
whereas Lemma (<ref>) is only stated for  and not for .
To get around this difficulty, we consider a larger set ⊃ such that
 is full dimensional and ϕ is area convex on .
Then we use Remark <ref> to obtain the final result, i.e., area convexity of ϕ.

	Let w = (x,u,y,z) and define 

	ϕ(w) _i = 1^p∑_j = 1^n P_ijγ_p_i(x_j, y_i) +∑_i=1^pγ_2(u, y_i)+∑_i = 1^c∑_j = 1^n C_ijγ_c_i(x_j, z_i) +∑_i = 1^cγ_2(u, z_i),
 where p_i = 2*P∞/P_i,:1 and c_i = 2*C∞/C_i,:1, then ϕ is area convex with respect to 1/3[   P -_p;  -C  _c ] on set :=_+,∞^n+1(1) ×Δ_p^+ ×Δ_c^+. In particular, it also implies 6√(3)ϕ is area convex with respect to 2√(3)[   P -_p;  -C  _c ] on set .





We also consider an δ-optimal solution oracle for ϕ as follows

	A δ-optimal solution oracle (OSO) for ϕ:→ takes input a and outputs w ∈ such that 
	
    a^Tw -ϕ(w) ≥sup_w∈ a^Tw -ϕ(w) -δ
 
 

Given Φ as a δ-OSO for function 6√(3)ϕ then consider the following algorithm:

For Algorithm <ref>, <cit.> shows the following:

	Let ϕ: →[-ρ, 0]. Suppose ϕ is area convex with respect to 2√(3)[   P -_p;  -C  _p ] on . Then for all t≥ 1 we have w_t/ t ∈ and, 
    sup_w∈w Jw_t/t≤δ + ρ/t

	In particular, in ρ/ iterations, Algorithm <ref> obtain (δ+)-solution of the saddle point problem (<ref>).

Theorem <ref> addresses the first part of the key three challenges.
Next, Lemma <ref>
shows an upper bound on the range of ϕ.

	Function ϕ:→ [-ρ,0] then ρ = O(P∞logp + C∞logc).


Finally, we need an efficient δ-OSO. Consider the following alternating minimization algorithm.

<cit.> shows the following convergence result.

	For δ > 0, Algorithm <ref> is a δ-OSO for ϕ which converges in O(log1/δ)
	iterations.

We show that for our chosen ϕ, we can compute the two argmax in each iteration of Algorithm <ref> analytically with computation time O((P)+(C)) and hence we obtain a δ-OSO running in O(((P)+(C)log1/δ) total work.
Parallelizing matrix-vector multiplications, eliminaters the dependence on (P) and (C),
at the cost of another log(N) term.

Each  in Algorithm <ref> can be computed as follows:

	x^k= min[]exp[]a/P^T y^k + C^T z^k -1, _n for all j ∈ [n].

	y^k = _Δ_p^+[]exp[]1/2(P∞ +1)(a^1-Px^k-1logx^k-1)

	z^k = _Δ_c^+[]exp[]1/2(C∞ +1)(a^2-Cx^k-1logx^k-1)

	In particular, we can compute x^k, y^k, z^k in O((P)+(C)) work and O(log N) parallel time.





§ ACKNOWLEDGEMENTS

We thank Richard Peng for many important pointers and discussions.

acm






§ PROOF OF AUXILIARY RESULTS


In this section, we include proofs of lemmas from
the main paper.
In some cases, the lemmas are direct restatements
of results from other papers, for which we provide appropriate pointers.



 §.§ Proof of Lemma <ref>

Suppose there exists x satisfying MPC.
Then, max_y ∈Δ^+_c,z ∈Δ^+_p  L(x, y,z) = 0.
On the other hand, if there does not exist any x satisfying MPC
then max_y ∈Δ^+_c, z ∈Δ^+_p  L(x,y,z) will indicate
the value of largest infeasibility among the packing-covering constraints
for the given x, and this value would be positive for all possible values of x
due to the infeasibility of this MPC instance.	



Consider an arbitrary strongly convex function d.
Assume WLOG that d(0)=0.
(otherwise, we can shift it accordingly).
We will show that max_x∈_∞^n(r)d(x) ≥nr^2/2 by induction on n for set _∞^n(r).
This suffices because _+,∞^n(1) is isomorphic to _∞^n(1/2). 
The claim holds for n = 1 by the definition of strong convexity.
Now, suppose it is true for n-1.
Then there exists x∈_∞^n-1(r) such that d(x) ≥(n-1)r^2/2.
Moving r units in the last coordinate from x in the direction of nonnegative slope, suppose we reach x∈_∞^n(r).
Then, due to strong convexity of d, we have d(x) ≥ d(x) + 1/2x-x∞^2 ≥(n-1)r^2/2 + r^2/2 = nr^2/2.



	Suppose we are given (x,y,z) such that sup_(x, y, z) ∈_+,∞^n ×Δ_p^+ ×Δ_c^+ L(x, y, z) - L(x, y, z) ≤.
	If there exists x which is feasible for MPC then choosing x = x then L(x, y, z) ≤ 0. Hence we have 
	
    sup_(y, z) ∈Δ_p^+ ×Δ_c^+L(x, y, z)    ≤
    ⇒Px - _p_+∞ + -Cx + _c_+∞   ≤,
 where implication follows by optimality over extended simplices Δ_p^+, Δ_c^+. So we obtain, if there exist a feasible solution for MPC then x is -approximate solution of MPC.

	On the other hand, suppose x is not an -approximate solution. Then 
	
    max   Px - _p_+∞, -Cx + _c_+∞ > 
    ⇒sup_(y, z) ∈Δ_p^+ ×Δ_c^+L(x, y, z) =    Px - _p_+∞ + -Cx + _c_+∞ >

	Let (y, z) ∈Δ_p^+ ×Δ_c^+ such that L(x, y, z) > then we have 
	
    sup_x∈_+,∞^n L(x, y, z) - L(x, y,z)     ≤
    ⇒ L(x, y, z) - inf_x∈_+,∞^n  L(x,y,z)    ≤
    ⇒inf_x∈_+,∞^n  L(x,y,z)  > 0

	Hence, if x is not -approximate solution of MPC then (y,z) satisfy y^T(Px -_p) + z^T(-Cx+_c) > 0 for all x∈_+,∞^n(1) implying that MPC is infeasible.




Let B = [  0 -1;  1  0 ]
and T := [  A -B;  B  A ].

Then A  B  iff T ≽ 0 iff all principle minors of T are nonnegative. Now, T ≽ 0 implies A ≽ 0. It is easy to verify that third principle minor is nonnegative iff (A) ≥ 1. So T ≽ 0 implies A must be invertible. Then, applying Schur complement lemma, we obtain that T ≽ 0 ⇔ A + BA^-1B ≽ 0. Now let A = [ a b; b d ] then A^-1 = 1/ad-b^2[  d -b; -b  a ]. It is easy to verify that A + BA^-1B = A(1-1/(A)). This implies T ≽ 0 ⇔ A ≽ 0 and (A) ≥ 1. Hence we conclude the proof.



This lemma appears exactly as Theorem 1.6 in <cit.>.
The proof follows from the same.




	
		
  * 
    A  B    ⇔[  A -B;  B  A ]≽ 0 
       ⇔ x^TAx + y^TAy + y^TBx -x^TBy  ≥ 0,   ∀ x, y
       ⇔ x^TAx + y^TAy  - y^TBx + x^TBy ≥ 0,  ∀ x, y
       ⇔[  A  B; -B  A ]≽ 0 ⇔ A  (-B)

		Here, the third equivalence follows after replacing y by -y. Hence we conclude the proof of part 1.
		
  * 
    A  B ⇔[  A -B;  B  A ]≽ 0 ⇒[  λ A -λ B;  λ B  λ A ]≽ 0 ⇔λ A ≽λ B

		
  * A  B implies [  A -B;  B  A ]≽ 0. Similarly C  D implies [  C -D;  D  C ]≽ 0. Hence 
    [    A+C -(B+D);  (B+D)  (A+C) ]≽ 0.
 So we obtain A+C  (B+D).
	




This lemma appears as Theorem 1.3 in <cit.>, and
the proof follows from the same.




We use equivalent characterization proved in Lemma <ref>. We need to show that d^2γ_β(a,b) ≽ 0 and (d^2γ_β(a,b)) ≥ 1 for all a ∈ (0,1] and b > 0. First of all, note that d^2γ_β is well-defined on this domain. In particular, we can write

    d^2γ_β(a,b) = [    β/b 1+loga; 1+loga    b/a ].

Note that a 2× 2 matrix is PSD if and only if its diagonal entries and determinant are nonnegative. Clearly diagonal entries of d^2γ_β(a,b) are nonnegative for the given values of β,a and b. Hence, in order to prove the lemma, it suffices to show that (d^2γ_β(a,b)) ≥ 1.

(d^2γ_β(a,b)) = β/a - (1+loga)^2 is only a function of a for any fixed value of β≥ 2. Moreover, it can be shown that (d^2γ_β) is a decreasing function of a on set (0,1]. Clearly, the minimum occurs at a =1. However, (d^2γ_β(1, b)) = β -1 ≥ 1 for all b > 0. Hence we have that (d^2γ_β(a,b)) ≥ 1 for all a ∈(0,1], b > 0 and β≥ 2.

Finally to see the claim that (d^2γ_β) is a decreasing function of a ∈ (0,1] for any β≥ 2, consider 

    d/da[](d^2γ_β(a,b))   = -β/a^2 - 2(1+ loga)/a
       ≤ -2(1+a(1+loga))/a^2 < 0
where the last inequality follows from the observation that 1 +a +aloga > 0 for all a∈ (0,1]. Hence we conclude the proof.




Note that γ_c_i, γ_p_i are twice differentiable in the (). So by Lemma <ref> part 2, it is sufficient to prove that d^2ϕ(w)  J for all w ∈().

By definition, we have γ_c_i≥ 2 for all i ∈ [c] and γ_p_i≥ 2 for all i ∈ [p]. Moreover x_j ∈ (0,1) and y_i > 0, z_i > 0 for any w = (x,u,y,z) ∈(). Then by Lemma <ref> and Proposition  <ref>, we have

    d^2ϕ(w)    = ∑_i = 1^p ∑_j = 1^n P_ij d^2 γ_p_i(x_j, y_i) + ∑_i=1^p d^2γ_2(u, y_i)+ ∑_i = 1^c ∑_j = 1^n C_ij d^2 γ_c_i(x_j, y_i) + ∑_i=1^c d^2γ_2(u, z_i) 
       (∑_i =1^p ∑_j = 1^n -P_ij e_j  ⊗ e_n+1+i + ∑_i=1^p e_n+1⊗ e_n+1+i
                + ∑_i = 1^c ∑_j = 1^n C_ij e_j ⊗ e_n+p+i + ∑_i=1^c (-1)e_n+1⊗ e_n+1+p+i)  ,

where e_k⊗ e_l = e_k e_l^T - e_le_k^T. Here we used P_ijd^2γ_p_i(x_j, y_i)  -P_ije_j⊗ e_n+1+i using Lemma <ref>, Proposition <ref> part 1, part 2 and C_ij d^2γ_c_i(x_j, y_i)  C_ij e_j ⊗ e_n+1+p+i using Lemma <ref>, Proposition <ref> part 2. Similar arguments can be made about terms inside the other two summations. Finally we used Proposition <ref> part 3 to obtain (<ref>). Note matrix in the last sum term is in fact J.

It is clear that since d^2 ϕ  J hence using Proposition <ref> part 2, we have d^2 6√(3)ϕ 6√(3) J. Then by Lemma <ref> part 2, we obtain 6√(3)ϕ is area convex with respect to 2√(3)[   P -_p;  -C  _c ] on set .

Note that the set of interest ⊂. Moreover,  is a convex subset. By Remark <ref>, one can see that 6√(3)ϕ is area convex with respect to 2√(3)[   P -_p;  -C  _c ] on set . Hence we conclude the proof.




Note that γ_β(a,b) ≤ 0 for any a ∈ [0,1], b ∈ [0,1], β≥ 0. Since P_ij≥ 0, C_kj≥ 0 for all possible values of i,j,k hence we clearly have ϕ(w) ≤ 0 for all w ∈. Now we prove that lower bound is not too small.

We have

    ∑_i =1^p∑_j = 1^n P_ijγ_p_i(x_j,y_i)    = ∑_i =1^p∑_j = 1^n P_ij*y_i x_jlogx_j + p_i y_i logy_i
       ≥ -∑_i =1^p∑_j = 1^n P_ij y_i 1/e + ∑_i=1^pp_iy_ilogy_i∑_j = 1^pP_ij
       = -∑_i =1^p∑_j = 1^n P_ij y_i 1/e + ∑_i=1^p2P∞y_ilogy_i
       ≥ -∑_i=1^p P∞/e y_i + ∑_i=1^p2P∞y_ilogy_i
       ≥ -P∞/e - 2P∞logp

Note that if w ∈ implies u  = 1. So 

    ∑_i =1^p γ_2(u,y_i) = ∑_i = 1^p 2y_ilog(y_i) ≥ -2 logp

Similarly, we have

    ∑_i =1^c∑_j = 1^n C_ijγ_c_i(x_j,z_i)    ≥ -C∞/e - 2C∞logc
    ∑_i =1^c γ_2(u,z_i)    ≥ -2 logc

Taking sum of all four terms, we conclude the proof.




Note that maximization with respect to u is trivial since u = 1 is a fixed variable. We first look at maximization with respect to x ∈_+,∞^n(1). Writing the first order necessary condition of Lagrange multipliers, we have

    a_j - ∑_i=1^p P_ijtγ_p_i(t,v)(t,v) = (x_j, y_i)- ∑_i=1^c C_ijtγ_c_i(t,v)(t,v) = (x_j, z_i) - λ_j    = 0 
    ⇒ a_j - []∑_i=1^p P_ij y_i + ∑_i=1^c C_ij z_i (1+logx_j) - λ_j    = 0.

Here λ_j is the Lagrange multiplier corresponding to the case that x_j = 1. By complimentary slackness, we have λ_j > 0 iff x_j = 1.

This implies x_j = min*exp*a_j/∑_i=1^p P_ij y_i + ∑_i=1^c C_ij z_i -1, 1 for all j ∈ [n].

Now we consider maximization with respect to y,z. Note that there are no cross-terms of y_i and z_i, i.e., γ_p_iy_i is independent of z variable and vice-versa. So we can optimize them separately. From first order necessary condition of Lagrange multipliers for y, we have

    a^1_i -∑_j = 1^n P_ijvγ_p_i(t,v)(t,v) = (x_j, y_i) -  vγ_2(t,v)(t,v) = (u, y_i) - λ   = 0
    ⇒ a^1_i - ∑_j =1^n P_ij (x_jlogx_j +p_i(1+logy_i)) - uloguu =1 - 2(1+logy_i) -λ   = 0
    ⇒ a^1_i - ∑_j =1^n P_ij x_jlogx_j  - 2(P∞+1)(1+logy_i) -λ   = 0

where last relation follows due to definition of p_i and λ is Lagrange multiplier corresponding to the constraint ∑_i=1^p y_i ≤ 1. By complimentary slackness, we have λ > 0 iff ∑_i=1^p y_i = 1.

Eliminating λ from above equations, we obtain y = _Δ_p^+*exp*1/2(P∞ +1)(a^1-Pxlogx).

Similarly, we obtain z = _Δ_c^+*exp*1/2(C∞ +1)(a^2-Cxlogx).

It is clear from the analytical expressions that for each iteration of Algorithm <ref>, we need O((P)+(C)) time. Hence total runtime of Algorithm <ref> is O(((P)+(C))log1/δ).




§ PROOF OF WIDTH REDUCTION FOR THE MPC PROBLEM
 

In Section <ref>, we made the assumption
that all entries

This assumption follows from the results in <cit.>.
We outline this proof in this section for completeness.

For the purpose of this proof, we introduce notation [k] := {1,…, k}.

Suppose we are given an instance of mixed packing covering of the form 

    Px ≤_p, Cx ≥_c, x≥_n.

Case 1: For each column P_:, i associated with variable x_i, let P_j_i,imax_j ∈ [p]P_ji > 0.
Then we consider the following updates to MPC in order to reduce diameter.

Suppose, without loss of generality, C_1,i = max_j ∈ [c] C_ji and C_ci = min_j ∈ [c]C_ji. If C_1i≤ P_j_i, i then we can update P_:,i = 1/P_j_i,iP_:, i, C_:,i = 1/P_j_i,iC_:, i and x_i = P_j_i,i x_i.
Then we observe that each element in P_:,i, C_:,i is at most 1.
Moreover, due to the packing constraint P_j_i,:x≤ 1,
we note that for any feasible x, P_j_i, ix_i ≤ 1. Finally, since P_j_i,i = 1, we have that x_i ≤ 1 lies in the support of constraint set. So we replaced the i-th column and corresponding i-th variable of the system by an equivalent system.

Similarly, if C_c,i≥ P_j_i,i then consider x^sol defined as

    x^sol_k:= 1/P_j_i,i   if  k = i
     0    otherwise.

Then x^sol is already a feasible solution of MPC.
So we may assume that C_ci < P_j_i,i < C_1i.
In this case, define r_i = C_1i/P_j_i,i and n_i = logr_i.
We make n_i copies of the column C:,i and denote by the tuple (i,l)
the columns of a new matrix C_:,(i,l) where l ∈ [n_i].
Similarly, we add n_i copies of variable x_i, denoted as x_(i,l).
We make similar changes to P_:,i.
Note that this system is equivalent to earlier system in the sense that any solution x_(i,l), l ∈ [n_i] can be converted into a solution of the earlier system
since x_i = ∑_l ∈ [n_i]x_(i,l).
However, this allows us to reduce the elements of C along with certain box constraints on x_i, which was our original goal.
For each j ∈ [c], l ∈[n_i], redefine 
    C_j,(i,l) = min{C_ji, 2^lP_j_i,i}
 
and for variable x_(i,l),
add the constraint

    x_i,l≤2/2^l P_j_i,i.
 

	MPC (<ref>) and the new system defined by matrices C, P and variable x are equivalent.


	For this proof, let us focus on i-th column and i-th variable.
	
	For any feasible solution x, consider x_i = ∑_l ∈ [n_i]x_i,l.
	This x_i does not violate any covering constraint since C_j, (i,l)≤ C_ji. 
	The packing constraints also follow because we have not made any changes to the elements
	corresponding to the packing constraints P_j, (i,l).
	
	For the other direction, the key fact to note
	is that any feasible x satisfies x_i ≤1/P_j_i,i
	due to packing constraint P_j_i,:x ≤ 1.
	Let l_i be the largest index such that 
	
    x_i ≤2/2^l_iP_j_i,i,

	and then let 
	
    x_(i,l) = 
    	x_i    if  l = l_i 
     0    otherwise.

	By construction, x_(i,l) satisfies the constraint in (<ref>)
	for all l ∈[n_i].
	Moreover, for constraint j, we must have C_j,:x≥ 1. 
	Note that if C_j,(i, l_i) = C_ji then there is nothing to prove.
	So we assume that C_ji > C_j,(i, l_i) = 2^l_i P_j_i, i.
	Then we must have that l_i < n_i in this case, by definition of n_i.
	This then gives x_(i,l_i) = x_i ≥1/2^l_i P_j_i,i
	by our choice of l_i being the largest possible.
	Then we know that C_j,(i,l_i) = 2^l_i P_j_i, i, and hence the
	j-th covering constraint is satisfied.
	
	Packing constraints are satisfied trivially since there is no change in elements of P_:,(i,l) for all l ∈[n_i]. Hence the claim follows.


Finally the proof follows by change of variables as x_(i,l) = 2^l-1 P_j_i,i and C_:,(i,l) = 1/2^l-1P_j_i,iC_:,(i,l).
Further, note that all elements of P_:, (i,l) are at most 1 for all l ∈ [n_i],
and all elements of C_:,(i,l) are at most 2 for all l ∈ [n_i] and x_i,l≤ 1 for all l ∈ [n_i].

Case 2:  Suppose P_j_i, i   = 0.. This implies that in variable x_i, this is a purely covering problem.
So we can increase x_i to satisfy the jth covering constraint such that
C_ji > 0 independent of the packing constraints and problem reduces to smaller packing covering problem in remaining variables and covering constraints j such that C_ji = 0. For this smaller packing covering problem, we can apply the method in Case 1 again. 






§ APPLICATION TO THE DENSEST SUBGRAPH PROBLEM
 

In this section, we apply the result in Theorem <ref>
to the densest subgraph problem.

We define the density of a graph G = ⟨ V,E ⟩
as |V|/|E| (half the average degree of G).
Hence, the densest subgraph
of G is induced on a subset of vertices U ⊆ V such
that

    U _S ⊆ V|E(S)||S|,

where E(S) denotes the set of edges in the subgraph of G
induced by S.

The following is a well-known LP formulation of the densest subgraph problem,
introduced in <cit.>,
which we denote using Primal(G).
The optimal objective value is known to be ρ_G^*.


    [          maximize           ∑_e ∈ E               y_e;        subject to                      y_e ≤ x_u,x_v,        ∀ e=uv ∈ E;                             ∑_v ∈ V          x_v ≤ 1,;                                     y_e ≥ 0, x_v ≥ 0,  ∀ e ∈ E, ∀ v ∈ V ]


We then construct the dual LP for the above problem.
Let f_e(u) be the dual variable associated with the first 2m
constraints of the form y_e ≤ x_u,
and let D be associated with the last constraint.
We get the following LP, which we denote by Dual(G),
and whose optimum is also ρ_G^*.


    [                minimize                                               D;              subject to                f_e(u) +             f_e(v) ≥ 1,              ∀ e=uv ∈ E;                                         ∑_e ∋ v             f_e(v) ≤ D,                 ∀ v ∈ V;                                                 f_e(u) ≥ 0, f_e(v) ≥ 0,              ∀ e=uv ∈ E ]


Parametrizing with respect to D,
this becomes a mixed packing covering LP.
The solution to the densest subgraph problem is
simply the smallest value of D for which the LP is feasible.
Since D can take at most O(|V||E|) ≤ O(|V|^3) values in total,
the densest subgraph problem can be reduced
to solving O(log |V|) instances of MPC,
where the number of nonzeros N in the matrix is O(|E|)
and the width w is simply the maximum degree in G.
This gives the following corollary.

 
Given a graph G = ⟨ V,E ⟩ with maximum degree Δ,
we can find the (1+)-approximation to the
maximum subgraph density of G, ρ_G^*,
in parallel time O(Δ^-1)
and total work O(Δ |E| ^-1).


The previous fastest algorithms for densest subgraph
do not depend on Δ - however,
their dependence on 1/ is quadratic <cit.>.
Corollary <ref> gives the fastest algorithm for this problem
in the high precision regime (< 1/Δ), since its dependence
on ^-1 is only linear.


