Multi-Dimensional Balanced Graph Partitioning via

Projected Gradient Descent∗

Dmitrii Avdiukhin Sergey Pupyrev Grigory Yaroslavtsev

Indiana University Facebook Indiana University
Bloomington, IN Menlo Park, CA Bloomington, IN
davdyukh@[Link] spupyrev@[Link] grigory@[Link]

ABSTRACT The process of splitting the input graph into these subsets,
Motivated by performance optimization of large-scale graph also known as graph partitioning, is essential for optimizing
processing systems that distribute the graph across multi- performance of such systems [3, 18, 20, 43].
ple machines, we consider the balanced graph partitioning Created partitions have a significant impact on the com-
problem. Compared to most of the previous work, we study munication between different workers and the resource usage
the multi-dimensional variant in which balance according to of individual workers. In order to maximize the process-
multiple weight functions is required. As we demonstrate by ing speed, the partitions should largely be independent to
experimental evaluation, such multi-dimensional balance is minimize communication. At the same time, computation
essential for achieving performance improvements for typical executed on each partition should take approximately the
distributed graph processing workloads. same amount of processing time, as the overall performance
We propose a new scalable technique for the multi- depends on the slowest worker. These constraints give rise to
dimensional balanced graph partitioning problem. It is based the Balanced Graph Partitioning problem whose goal
on applying randomized projected gradient descent to a non- is to divide the vertices of a graph into a given number of
convex continuous relaxation of the objective. We show how (approximately) equal size components while minimizing the
to implement the new algorithm efficiently in both theory resulting edge cut. Balanced Graph Partitioning is a
and practice utilizing various approaches for the projection classic and thoroughly studied problem from both theoretical
step. Experiments with large-scale graphs containing up to and practical points of view [9, 12]. In the context of dis-
hundreds of billions of edges indicate that our algorithm has tributed graph processing, the problem is typically studied
superior performance compared to the state of the art. in two variants.
In the vertex partitioning model, each worker machine is
PVLDB Reference Format: assigned an equal number of vertices with the goal of mini-
Dmitrii Avdiukhin, Sergey Pupyrev, Grigory Yaroslavtsev. Multi- mizing the number of cross-machine edges. Since messages
Dimensional Balanced Graph Partitioning via Projected Gradient are usually sent between adjacent vertices, storing tightly
Descent. PVLDB, 12(8): 906 - 919, 2019.
DOI: [Link] connected subgraphs on the same worker can reduce com-
munication and hence running times of jobs. It has however
been observed that this strategy does not lead to equally
1. INTRODUCTION loaded partitions for real-world graphs with power law de-
Distributed graph processing systems have been widely gree distribution [18]. Graph partitioning algorithms tend to
adopted in recent years to enable analysis and knowledge ex- collocate high-degree vertices and corresponding partitions
traction from large-scale graphs. Systems such as Giraph [6], take much longer to process, resulting in longer execution
GraphX [19], GraphLab [30], and PowerGraph [18] allow time overall.
users to use a vertex-centric model for applications which The edge partitioning model has been suggested to alleviate
can be executed on a cluster of worker nodes. In this setting, the above imbalance problem [18, 29]. In this model the goal
each worker node operates on a subset of the input graph is to partition the graph so that the number of edges in every
and communicates with other workers by sending messages. component is the same, while the number of incident edges
across different components is minimized. Good partitions
∗Research supported by NSF award CCF-1657477 and Face- according to this model typically result in better balance
book Faculty Research Award. Full version: [Link] across workers and reduced computation time in comparison
org/abs/1902.03522. to the trivial hash-based assignment of vertices to worker
machines. However, edge-based graph partitioning can still
result in performance regressions [3, 40].
This work is licensed under the Creative Commons Attribution- To analyze the source of regressions, we performed a simple
NonCommercial-NoDerivatives 4.0 International License. To view a copy experiment of running a Page Rank algorithm implemented
of this license, visit [Link] For
any use beyond those covered by this license, obtain permission by emailing
on top of Giraph utilizing various graph partitioning meth-
info@[Link]. Copyright is held by the owner/author(s). Publication rights ods. Figure 1 illustrates the histograms of running times for
licensed to the VLDB Endowment. individual workers processing a graph with 800M vertices
Proceedings of the VLDB Endowment, Vol. 12, No. 8 and 80B edges. As discussed above, partitions according to
ISSN 2150-8097. the vertex partitioning model suffer from unequal distribu-
DOI: [Link]
 40

iteration time, sec

30

20

10

0
hash vertex partitioning edge partitioning vertex-edge partitioning
6.25% 72.72% 68.84% 59.23%

Figure 1: The running time of an iteration of Page Rank on a Giraph cluster of 16 worker machines using various graph
partitioning strategies. The numbers indicate the average percentage of local (uncut) edges per worker, which is proportional to
the number of local messages for the distributed graph processing workload. Vertex-edge partitioning achieves approximately
12% iteration time improvement compared to hash.

tion of edges across workers. A single overloaded partition as Multi-Dimensional Balanced Graph Partitioning
contains 1.92x more edges than an average one, which re- (MDBGP).
sults in 1.5x longer execution time. We also observe a high The simplest example of MDBGP is the classic balanced
correlation (ρ = 0.79) between the number of edges assigned graph partitioning problem, which is equivalent ot the ver-
to a partition and the corresponding processing time in this tex partitioning strategy described above and can be ex-
experiment. Partitioning according to the edge partitioning pressed using a single weight function w(1) (v) = 1. Since
model yields a 1.08x running time improvement over the w(1) (V ) = |V | this requires that we maximize edge locality
baseline, though there is still a noticeable imbalance between while ensuring that |Vi | ≈ |Vk | . Using two weight functions,
the fastest and the slowest worker machines. This can be
w(1) (v) = 1 and w(2) (v) = deg(v), corresponds to requiring
explained by uneven distribution of vertices among workers.
balance on the number of vertices and edges in the parts of
Machines with more vertices have higher operational over-
the partition and hence corresponds to the vertex-edge par-
head such as serialization of sent messages whose number is
titioning approach mentioned above. Indeed, w(2) (V ) = 2|E|
proportional to the number of vertices on a worker. Here we
and hence in additionP to balance on the number of vertices
observe an 1.33x imbalance in the number of vertices and a
this requires that v∈Vi deg(v) ≈ 2|E| k
. However, the model
moderate correlation (ρ = 0.62) between the running time
is not restricted to vertex- and edge-balance (as in the afore-
and the vertex count on the workers.
mentioned vertex-edge partitioning) but can take arbitrary
In order to mitigate the issues described above we in-
user-specified weights. In particular, when partitioning the
troduce a new strategy, vertex-edge partitioning, which is
vertices of the graph between the workers for load balancing,
designed to balance the number of vertices and edges across
various weights modeling expected vertex activity can be
workers simultaneously. As shown in Figure 1, this is done at
used (historical data on individual vertex load, proxy values
a cost of a lower edge locality (percentage of edges with both
for the load such as PageRank, etc).
endpoints on the same machine), and thus, higher communi-
While a large body of work exists offering practical
cation volume. The resulting assignment results in a 1.12x
solutions for the one-dimensional version of the prob-
speedup over the hash-based model. Motivated by the above
lem [7, 12, 13, 14, 22, 23, 33, 41, 42], as well as on theoretical
experiment and a number of earlier studies [3, 20, 40, 43],
foundations of graph partitioning [4, 26, 32], literature on
we formalize a new model for graph partitioning, which is
principled and scalable approaches for the multi-dimensional
suitable for real-world distributed graph processing systems.
case is quite sparse [24, 35, 36, 37]. In particular, if the weight
We now formally describe the model in the most general
functions are unrelated to each other, one can easily con-
setting which allows one to require balance according to d
struct examples when no feasible solution exists that satisfies
different unrelated weight functions. Let G(V, E) be a graph
all balance constraints even for two weight functions. How-
with d vertex weight functions w(1) , . . . , w(d) : V → R+ ,
ever, it is empirically observed that instances coming from
each assigning aP positive weight to every vertex in the graph.
applications often allow balanced solutions for several weight
Let w(j) (V ) = v∈V w(j) (v) be the sum of weights of all
functions of interest simultaneously. For classical local search
vertices in the graph according to the j-th weight function.
based algorithms such as [25], handling of multiple unrelated
Given an integer k and a parameter ε > 0, the goal is to
weight functions is challenging since imposing one balance
partition the vertex set V into k sets V1 , . . . , Vk such that,
constraint might violate another; hence, finding a good local
for each weight function w(j) and each part Vi , the sum
move becomes computationally hard. We overcome this dif-
of weights in Vi is approximately the same and close to
(j) ficulty by using a continuous relaxation of the problem. In
(j)
(v) = (1 ± ε) w k(V ) . We
P
the average, that is, v∈Vi w order to obtain an integral solution, we apply randomized
call such partitions ε-balanced. Finally, among all such - rounding which preserves balance with high probability.
balanced partitions the goal is to find one that maximizes
the number of edges whose both endpoints are contained
within some component of the partition and hence min- 1.1 Our Contributions
imizes the size of the cut. This problem is referred to We present a scalable algorithmic framework for the prob-
lem of balanced partitioning of large graphs according to
multiple user-specified weight functions while maximizing the (METIS [23], Fennel [41]), branch-and-bound [13], label prop-
number of edges inside the resulting components. Our frame- agation and local search (Balanced Label Propagation [42],
work consists of applying the projected gradient descent on a Social Hash Partitioner [22], Spinner [33]), as well as hybrid
standard relaxation with a suitably chosen projection method. approaches (linear embedding method combined with various
The relaxation is to maximize a non-convex quadratic func- optimizations [7]). Due to the combinatorial nature of these
tion f (x) = 12 xT Ax for x ∈ Rn , where A is the adjacency algorithms, their generalizations to the multi-dimensional
matrix, subject to a constraint x ∈ K for a certain convex case appear to be non-straightforward without substantial
body K defined by the weight functions. Section 2 provides losses in performance, while our continuous relaxation han-
the exact description of the relaxation. Note that the gra- dles multiple balance constraints uniformly. Compared to
dient descent step only uses a matrix-vector multiplication the one-dimensional version, existing literature on the multi-
since ∇f = Ax. Thus, the algorithm allows a straightforward dimensional version is rather sparse [24, 35, 36, 37] and the
distributed implementation. main publicly available tool for the problem is currently
While applying projected gradient descent to solve non- METIS [24, 37].
convex optimization problems subject to convex constraints Vast literature exists on optimization of non-convex func-
is a well-studied approach in non-linear optimization (Section tions and the interest in this topic lately has been particularly
2.3, [8]) and machine learning (Section 6.6, [21]), one has to high. However, in the constrained case when the optimization
overcome two technical challenges to make it applicable to the has to be performed over a convex body, fairly little is known;
multi-dimensional graph partitioning problem: 1) projection see classic optimization literature [8,11,44]. Recent results on
step is computationally expensive, 2) existence of points with the non-convex optimization problem subject to convex con-
small gradient (saddle points) slows down convergence. straints and its special cases include [5, 16, 17, 21, 39]. Closest
We show how to address the first challenge by designing to our work in terms of techniques is [27] who use projected
efficient projection step algorithms tailored to the standard gradient method to solve convex programs involving the max-
relaxation of MDBGP. While convergence to the projection norm and show how to solve large semidefinite programming
point can be achieved using various alternating projections relaxations of Max-Cut. Their results are quite different
methods [15], for d ≤ 2 we give one-shot exact solutions with from ours as we consider a balanced version of graph parti-
almost linear running time. tioning and expect our algorithms to be scalable; the largest
instances handled by [27] have |V | = 20K and |E| = 40K.
Theorem 1.1. Running time of the projected gradient Since we require that our algorithms scale to graphs with
descent step is O(|E| + |V | logd−1 |V |) for d ≤ 2 and scales as billions of edges, using existing general purpose software for
O(|E|/m+|V | logd−1 |V |) when distributed over m machines. constrained quadratic programming is infeasible.

In order to address the second challenge, we use small per- 2. PROJECTED GRADIENT DESCENT
turbations to get out of saddle points, where the perturbation For an integer t ≥ 1, we use notation [t] to denote the
vectors are sampled from a scaled n-dimensional Gaussian set {1, . . . , t}. The weighted d-dimensional balanced graph
distribution. We refer to the resulting algorithm as Gradi- partitioning problem is defined by a collection of d weight
ent Descent (GD), see Algorithm 1. Convergence analysis of functions w(1) , . . . , w(d) , where w(j) : V → R+ . For a set
GD remains an open problem. While noisy gradient descent S ⊆ V we use notation w(j) (S) ≡ v∈S wv .
P (j)
is known to have fast convergence to a local optimum for
non-convex optimization subject to equality constraints, if Definition 2.1 (MDBGP). For a graph G(V, E), the
inequality constraints are allowed, convergence analysis is Multi-Dimensional ε-Balanced Graph k-Partitioning
unknown [16]. problem is to partition the vertex set V into k sets
Our experiments show that GD scales to graphs with up to V1 , . . . , Vk such that for each j ∈ [d], it holds that
(j)
several billions of vertices and several hundreds of billions w(j) (Vi ) = (1 ± ε) w k(V ) for all i ∈ [k]. Among all such
of edges. We perform a comparison of various graph parti- partitions the goal is to find one that maximizes the number
tioning strategies for optimizing several real-world Giraph of edges whose both endpoints are within some part.
workloads. The results demonstrate that multi-dimensional
balancing is a suitable objective for achieving performance In this paper we focus on the 2-partitioning problem; for
improvements, leading to speedups in the order of 10% − 30% the general variant of k-partitioning, we apply the algorithm
over the state-of-the-art one-dimensional strategies. Com- recursively. For k = 2 MDBGP is equivalent to the following
pared to existing scalable graph partitioners, such as Social integer quadratic program:
Hash Partitioner [22], Spinner [33], and Balanced Label Prop-
agation [34, 42], the algorithm is conceptually simple and
obtains close-to-perfect balance across multiple dimensions. 1 X
Maximize: (xi1 xi2 + 1)
2
(i1 ,i2 )∈E
1.2 Previous Work  n  n
X 
While one-dimensional balanced graph partitioning has (j) 
X (j)
Subject to: wi xi  ≤ ε wi ∀j ∈ [d]


been studied extensively and a number of tools exist [7,13,14, 
i=1

i=1
22,23,33,41,42] (see also surveys by Bichot and Siarry [9] and
xi ∈ {−1, 1} ∀i ∈ V
by Buluç et al. [12]), to the best of our knowledge none of the
practical algorithms for this problem have been previously
based on running gradient descent on a continuous relaxation. The interpretation of xi variables is that if xi = 1 then
Existing approaches are inherently discrete and are based on i ∈ V1 and if xi = −1 then i ∈ V2 . The objective is then
combinations of various discrete algorithms: greedy heuristics the same as in MDBGP and counts the number of edges
whose both endpoints are contained in some part of the Algorithm 1: GD (d-Dimensional Balanced Graph 2-
partition. Indeed, an edge (i1 , i2 ) makes a contribution of Partitioning via Randomized Projected Gradient Descent)
1 to the objective when xi1 = xi2 (and hence xi1 xi2 = 1)
input : Graph G(V, E), ε ∈ [0, 1], weight functions
and 0, otherwise (since xi1 xi2 = −1). The constraints are
w1 , . . . , w d : V → R +
equivalent to −εw(j) (V ) ≤ w(j) (V1 ) − w(j) (V2 ) ≤ εw(j) (V ). I−1
parameters: I, {ηt }t=0 , {γt }I−1
t=0
Adding or subtracting w(j) (V ) to both sides and dividing by
(j) output : ε-balanced partition w.r.t w(1) , . . . , w(d) of V
2, we have w(j) (Vi ) = (1 ± ε) w 2(V ) as required in MDBGP. into (V1 , V2 )
After dropping the additive term, the objective can be 1 K = B∞ ∩ dj=1 Sεj
T
expressed as f (x) = 12 xT Ax and has gradient ∇f (x) = Ax
2 x(0) = 0;
and Hessian ∇2 (f ) = A. Finally, we use a continuous relax- 3 for t = 0 to I − 1 do
ation of the above problem where we replace the integrality
4 z(t) = x(t) + Nn (0, ηt ); // Noise addition step
constraints with xi ∈ [−1, 1] for all i ∈ V . A solution to this
continuous relaxation can be converted into an integral solu- 5 y(t+1) = (I + γt A)z(t) ; // Gradient descent step
tion using randomized rounding. Using independent random 6 x(t+1) = argmin ky(t+1) − xk2 , // Projection step
x∈K
variables Xi for each vertex such that Pr[Xi = 1] = 1+x 2
i

1−xi
and Pr[Xi = −1] = 2 the expected value of the objec- 7 V1 = V2 = ∅; // Randomized rounding
tive on the rounded solution (X1 , . . . , X|V | ) is the same as 8 for each i ∈ V do
(t)
x +1
on the initial fractional solution (x1 , . . . , x|V | ) while all bal- 9 With probability i 2 , let V1 = V1 ∪ {i};
ance constraints are still approximately preserved with high otherwise, V2 = V2 ∪ {i};
probability by concentration bounds.

2.1 Overview
We propose the following algorithm for the multi- In the projection step of GD (Line 6) we need to find
argminx∈K ky(t+1) − xk2 , where K = B∞ ∩ dj=1 Sεj . Denot-
T
dimensional balanced graph partitioning problem based on
the continuous relaxation described above. The algorithm ing y(t+1) as y we formulate this step as an optimization
is referred to as Gradient Descent (GD), see Algorithm 1. It problem:
(t)
computes a sequence of vectors { x(t) }, where xi ∈ [−1; 1]
(0)
for all i ∈ V and t. Here x is initialized with zero vec-
tor, and x(t+1) is computed by applying projected gradient Minimize: f (x) = kx − yk22
descent iteration to x(t) . Each iteration consists of three Subject to: gi = x2i − 1 ≤ 0 ∀i ∈ [n]
steps. n
Step 1: Adding noise. We add Gaussian noise to x(t) ob- (j)
h+ =
X (j)
wi xi − ε ≤ 0 ∀j ∈ [d]
taining a noisy vector z(t) . The noise is drawn from the n- i=1
dimensional Gaussian distribution Nn (0, ηt ) with zero mean n
and variance ηt in each coordinate. The addition of noise to (j)
X (j)
h− = − wi xi − ε ≤ 0 ∀j ∈ [d]
x(t) allows to escape from saddle points, e.g. x(0) = 0. i=1
Step 2: Gradient descent. We obtain y(t) from the noisy
vector z(t) via a gradient descent step with step size γt . Note
that the gradient at z(t) is given as Az(t) hence this step can The optimum solution to the optimization problem has to
be expressed as y(t) = (I + γt A)z(t) . satisfy KKT conditions:
Step 3: Projection. The resulting Tdvector y(t) is then pro-
j
jected on the feasible space B∞ ∩ j=1 Sε , where: Stationarity:
n d n
B∞ = { x ∈ Rn | ∀i : xi ∈ [−1, 1] }
X X (j) (j)
X (j)
y−x= µi xi ei + (µ+ − µ− ) wi ei
n
X n
X i=1 j=1 i=1
(j) (j)
Sεj = { x ∈ Rn | | wi xi | ≤ ε wi } for j ∈ [d],
Complementary slackness 1:
i=1 i=1
µi (x2i − 1) = 0, ∀i ∈ [n]
that is, B∞ satisfies that kxk∞ ≤ 1 and Sεj corresponds to
the constraints imposed by the balance of weights according Complementary slackness 2:
to the j-th weight function. n
!
The final solution is obtained by rounding last x(t) : each (j)
X (j)
(t)
µ+ wi xi − ε = 0, ∀j ∈ [d]
x +1
vertex i is assigned to part V1 with probability i 2 . Note i=1
that this ensures that the expected number of edges whose n
!
(j)
X (j)
endpoints belong to the same part after this rounding is µ− wi xi +ε = 0, ∀j ∈ [d]
(t) (t)
given as 21 (i1 ,i2 )∈E (xi1 xi2 + 1).
P i=1

The algorithm uses parameters ηt , γt , and I, where t is the

(j) (j)
iteration index. Here ηt controls the magnitude of noise, γt Here µi , µ+ , µi ≥ 0 are the dual variables and ei is the
is the step size, and I is the number of iterations. We discuss i-th standard unit vector. It is a standard fact (see [11],
the selection of parameters in the experimental Section 4. Chapter 5.5.3) that for convex optimization subject to linear
constraints Stationarity, Complementary slackness and Pri-
2.2 Projection mal/Dual feasibility are necessary and sufficient conditions
for the optimum solution. Thus we just focus on satisfying y
these conditions
P below.
(j) (j) (j)
Let γi = dj=1 (µ+ − µ− )wi . Then by Stationarity for
each i we have yi − xi = µi xi + γi . Consider the following x
three cases: y − λ∗ w w
Case 1. (yi > 1 + γi ). If µi = 0 then by Stationarity
xi = yi − γi > 1 which violates primal feasibility conditions.
Therefore µi > 0 and x2i = 1 by Complementary slackness 1.
Among the two roots xi = 1 and xi = −1 the second root hw, xi = ε
can be ruled out and hence xi = 1. Indeed, if xi = −1 then
by Stationarity yi + 1 = −µi + γi which contradicts µi > 0
and yi > 1 + γi . hw, xi = −ε
Case 2. (yi < −1 + γi ). This case is symmetric to the
previous one and thus xi = −1 in this case.
Figure 2: One-dimensional projection. First, the initial
Case 3. yi ∈ [−1+γi , 1+γi ]. First we show that µi = 0. In-
point y is moved by vector −λ∗ w, which is an orthogonal
deed, assume that µi > 0. Then xi = ±1 by Complementary
vector to planes, corresponding to balance constraints. Then
slackness 1. Both cases lead to contradiction:
the resulting point is projected on the cube.
1. (xi = 1). By Stationarity yi − 1 = µi + γi which
contradicts with yi ≤ 1 + γi and µi > 0.
2. (xi = −1). Similarly to the above by Stationarity we for all j ∈ S− . In the analysis below we assume that
have yi + 1 = −µi + γi which is a contradiction with d = |S+ ∪ S− | corresponds to the “effective dimension” of the
yi ≥ −1 + γi and µi > 0. problem.
Therefore in this case we have µi = 0 and hence by Station-
arity xi = yi − γi . 2.3 Exact Projection Algorithms
(j) (j)
Let λj = µ+ − µ− and assume that these values are
known to the algorithm. For z ∈ R we use notation [z] = Projection for d = 1. As a warm up, we first show how
min(1, max(−1, z)) for the truncated linear function. Using to perform exact projection for d = 1 in O(n log n) time,
the analysis above the projection step is simply xi = [yi − proving Theorem 1.1 for d = 1. This can be further improved
Pd (j)
j=1 λj wi ]. It remains to show how to find { λj }. to O(n) using a more careful approach [31]. However, to the
Note that from Complementary slackness 2 it follows that best of our knowledge, no fast algorithm is known for d > 1
(j) (j)
either µ+ = 0 or µ− = 0 since both of these values being which is the main focus of our work. Dropping the second
(1)
positive leads to a contradiction. This leads to three cases: 1) index to simplify presentation (that is, wi = wi ) and using
(j) (j) (j) (j) (j) (j)
µ+ = 0, µ− > 0, 2) µ− = 0, µ+ > 0 and 3) µ+ = µ− = 0 the fact that xi = [yi − λwi ] we have:
which correspond to the three possibilities for sign(λj ). For
each of the d dimensions we can try all three choices. For X X X
a fixed guess of the signs let S+ = {j : λj > 0}, S0 = wi xi = wi + −wi +
{j : λj = 0} and S− = {j : λj < 0}. Assuming a correct i i : yi ≥1+λwi i : yi ≤−1+λwi
guess of sign(λj ) for each of the dimensions the optimization X
problem above reduces to the following: + wi (yi − λwi ).
i : yi ∈(−1+λwi ,1+λwi )
Proposition 2.1. For the correct guess of sign(λj ) for
all j ∈ [d] it suffices to find the optimum of the above opti- We introduce notation hi (λ) where each hi is the following
mization problem without the constraints for j ∈ S0 . This piecewise linear function:
optimum is unique.

The proof is given in the full version. Using Proposition 2.1 wi
 if λ < (yi − 1)/wi
and trying all guesses for sign(λj ) we can reduce the pro- hi (λ) = wi (yi − λwi ) if λ ∈ [(yi − 1)/wi , (yi + 1)/wi ]
jection step to 3d instances of the following optimization 
−w if λ > (yi + 1)/wi
i
problem:

Thus n
P Pn
i=1 wi xi = Pi=1 hi (λ) and the problem reduces to
Minimize: f (x) = ky − xk22 finding λ such that i hi (λ∗ ) = ±ε where the sign depends
∗

on whether our dimension is in S+ or S− . Since wiP ≥ 0 for

Subject to: gi = x2i − 1 ≤ 0 ∀i ∈ [n] all i each hi is monotone in λ and so the function i hi is
Xn
(j)
X n
(j)
a monotone piecewise linear function. The value of λ∗ can
wi xi = ε, ∀j ∈ S+ ; wi xi = −ε, ∀j ∈ S− be found in O(log n) iterations of binary search where each
i=1 i=1 iteration requires O(n) time to evaluate the sum. This gives
the overall running time of O(n log n). See Figure 2 for an
which can be done by finding numbers λj > 0 for illustration.
jP∈ S+ and λj < 0 for j ∈ S− and setting xi = [yi −
j∈S+ ∪S− λj wij ]. The choice of λj ’s has to satisfy the con- Projection for d = 2. For d = 2 we need to find (λ1 , λ2 ) such
(j) (j)
straints n
P Pn Pn
i=1 wij xi = ε for all j ∈ S+ and i=1 wij xi = −ε that i=1 hi (λ1 , λ2 ) = ±ε for j = 1, 2, where hi (λ1 , λ2 )
 y most expensive, in our experiments we mostly use the al-
ternating projections method. Moreover, since in practice
each iteration of alternating projection is computationally
y − λ1 w(1)
y − λ2 w(2) expensive, in the intermediate iterations we project on each
plane and the cube only once, while in the last iterations
x we run the alternating projections method until convergence.
(1) (2)
We refer to this choice as “one-shot” alternating projection
y − λ1 w − λ2 w below.

Figure 3: Two-dimensional projection. Initial point y is 3.2 Adaptive Step Size

moved by vector λ1 w(1) + λ2 w(2) and then projected on the Recall that Algorithm 1 has the following parameters:
cube. Gaussian noise variances for each step { ηt } and step size
parameters { γt }. Due to the spectral properties of the ad-
Table 1: Theoretical properties of projection methods. jacency matrix in our experiments, the algorithm does not
encounter any saddle points other than the initial point x = 0.
d Output Time required Therefore it suffices to add Gaussian noise only at the first
Alternating any x∈K Until convergence iteration (that is, ηt = 0 for t 6= 0).
Dykstra’s any projection Until convergence The simplest choice of the step size parameters { γt } is
Exact (ours) d≤2 projection O(n logd−1 n) constant, but it gives suboptimal results in our experiments.
Carefully chosen step size parameters for different iterations
not only gives better performance but can also be used to
ensure that convergence can be reached in a fixed number of
is defined below. steps. In section 4.3 we discuss how to choose the step size
 (j) to achieve good performance for a wide range of graphs.
wi
 if σi < yi − 1 The choice of step size parameters is complicated by the
(j)
hi (λ1 , λ2 ) = wi(j) (yi − σi ) if σi ∈ [yi − 1, yi + 1] projection step. The change in the objective function and the


−wi
(j)
if σi > yi + 1 progress towards an integral solution can both be related to
the progress in Euclidean distance kxt − xt+1 k between the
(1) (2)
where σi = λ1 wi + λ2 wi . The projection process is shown iterations, which can be greatly reduced by the projection.
in Figure 3. In Appendix A, we prove Theorem 1.1 for d = 2 Another important implementation detail is our handling
showing that nested binary search can be used to solve this of vertices which are close to integral. When the number of
problem in O(n log n) time. such vertices becomes large, the progress of the algorithm
can slow down. This is due to the fact that while the gradient
vector is still large, all of its large components correspond to
3. IMPLEMENTATION already integral vertices and point to the outside of the feasi-
ble region. These components dominate in the computation
3.1 Projection algorithms of the projection step, which leads to slow convergence. In
We considered the following three methods for the projec- order to avoid this issue, we “fix” such vertices so that they
tion step (Algorithm 1, Line 6). Their theoretical properties become integral and no longer participate in the gradient
are summarized in Table 1. update and the projection step. As we show in Section 4.3,
• Alternating projections: A standard approach for projec- this yields noticeable improvements in the quality of the
tion on the intersection of convex sets is the alternating results.
projections method (see [10]). It is easy to implement
projections on B∞ and ∩dj=1 Sεj separately. Since both
are convex bodies by alternating projections on each of 4. EXPERIMENTS
them one can guarantee convergence to a point in the We design our experiments to understand how well the new
intersection, but there is no guarantee that this point partitioning algorithm behaves on real-world datasets and
will be the actual projection. In practice, we are able how it affects the performance of distributed graph process-
to achieve slightly better balance by modifying this ing. We are primarily interested in two quality measures of
approach slightly and projecting on S0j instead of Sεj . partitions: cut size (number of non-local edges) and maximum
This still ensures that we get a point in the intersection imbalance over all dimensions (maximum relative deviation
in the end. of the weight in a part of the partition from the average). As
• Dykstra’s projection: We also consider Dykstra’s pro- pointed out in Section 1, we are not aware of an alternative
jection algorithm [1, 15]. This is a modification of the scalable approach for solving the multi-dimensional balanced
alternating projections method which is guaranteed to partitioning. However, some of the existing techniques for
converge to the projection. one-dimensional partitioning can be adapted for the multi-
• Exact projection for d ≤ 2: This is the algorithm pre- dimensional case. Next we discuss several such techniques
sented in Section 2.2. In our experiments Dykstra’s al- that we evaluated.
gorithm and exact projection give similar results, since Hash is the simplest partitioning strategy that assigns
they find approximately the same projection point. vertices to worker machines by hashing the vertex identi-
In Section 4.3 we study how quality of results produced by fiers. Hashing is stateless, extremely fast in practice, and
GD depends on choice of one of the projection methods above. requires no preprocessing of the graph, which made it the
Since the exact projection algorithm is computationally the default strategy in Giraph. The main disadvantage is that
the majority of sent messages are non-local and may results Spinner BLP SHP
0.5

0.41
0.39

0.81

0.51
in significant communication.

1.66
vertex imbalance
0.4
Spinner is a graph partitioning algorithm that can be

0.25
0.23
applied to process large-scale graphs in a distributed environ- 0.3
ment [33]. The algorithm is based on the label propagation 0.2

0.11
0.09

0.09
0.07
technique in which vertices exchange their labels trying to

0.03
0.1

0.02

0.02

0.02

0.02
pick the most frequent label among its neighbors. This pro-

0
0.0
cess guarantees a high number of adjacent vertices having the
0.10
same label, which are then assigned to the same worker. Spin-

edge imbalance
0.08
ner does not enforce a strict balance across partitions but

0.05

0.05

0.05

0.05

0.05
0.06

0.04
integrates score functions that penalize imbalanced solutions.
0.04

0.02

0.02
BLP is another approach based on the balanced label propa-

0.01
0.02
gation and combines the ideas of Ugander and Backstrom [42]

0
0

0
0

0
and Meyerhenke et al. [34]. On the first step, the method 0.00
k=2 k=8 k=2 k=8 k=2 k=8
creates a size-constrained clustering of the input graph us- LiveJournal Twitter Friendster
ing significantly more clusters than the number of available
machines, k. In our implementation, we construct c × k clus-
Figure 4: Vertex and edge imbalance ( max i w(Vi )
avgi w(Vi )
− 1) of
ters for c = 1024 and forbid a cluster to contain more than
|V | |E|
vertices and c×k edges. On the second step, we ran- the solutions created by different algorithms on the three
c×k
domly merge the clusters into k parts, which results in the public networks with k ∈ {2, 8}. Lower values correspond to
multi-dimensional balance even if the original clusters have more balanced partitions. Hash and GD yield near-balanced
different sizes. solutions for the instances.
SHP is a distributed graph partitioner [22, 38] that is based
on a classical local search heuristic [25]. Although SHP does
not provide balancing on multiple dimensions, it supports component affects the job performance. Figure 4 illustrates
a mode with several dimensions whose final balance is not the resulting vertex and edge imbalance of the solutions on
guaranteed. The algorithm works by balancing on a new the public networks for three algorithms, Spinner, BLP, and
dimension, which is a combination of the specified dimensions. SHP,
 using k =2 and k = 8. The imbalance is defined as
We configure SHP to find solutions having the same number maxi w(Vi )
− 1 , where the maximum and the average are
avgi w(Vi )
of edges (with a higher coefficient in the combination) and
taken over the total weight of all k constructed parts. Note
the same number of vertices (with a lower coefficient) in
that this definition differs from the usual notion of imbalance
every part.
in the sense that it considers only positive imbalance, since we
We implemented the algorithms and extensively experi-
are only interested in the most overloaded machines. We do
mented with the Giraph framework, which is used as the
not include the results for Hash and GD, as the corresponding
primary tool for large-scale graph analytics at Facebook [2,6].
values are below 0.01 for the instances.
Although the evaluation is performed with the single dis-
We observe that two algorithms, Spinner and SHP, are not
tributed graph processing system, we believe that our main
suitable for the multi-dimensional variant of the problem. For
conclusions are valid for other frameworks relying on the
dense graphs with a highly skewed degree distribution (as
vertex-centric programming model. For our experiments,
in Twitter), the algorithms cannot simultaneously provide
we use four large social networks that are publicly avail-
balance on the two dimensions. With the default setting,
able [28]. LiveJournal, Orkut, Twitter, and Friendster
these two algorithms generate solutions in which some of the
are undirected graphs containing 4.8, 3.1, 41, and 65 mil-
parts contain 1.5 − 2x more vertices than the average one. We
lion of vertices and 0.04, 0.12, 1.2, and 1.8 billion of edges,
tried to modify the techniques by adjusting relative weights of
respectively. In addition, we experiment with several large
their penalty functions for vertex and degree counts. However,
subgraphs of the Facebook friendship graph that serve to
we were not able to design universal penalty weights that work
demonstrate scalability of our approach and its performance
for all instances. A similar behavior regarding the resulting
on real-world data. We denote the graphs by FB-X, where X
balance is observed for larger graphs, FB-3B, FB-80B, and
indicates the (approximate) number of billions of edges; this
FB-400B. In contrast, Hash, GD, and BLP produced nearly-
data is anonymized before processing.
balanced (that is, having ε ≤ 0.05 both for vertex and edge
Next we analyze the quality of the solutions produced by
counts) solutions for all the instances. With this in mind, we
the algorithms on our dataset (Section 4.1) and evaluate var-
exclude Spinner and SHP from further experiments.
ious graph partitioning strategies for speeding up distributed
Next we compare the quality of our algorithm as measured
graph processing for real-world workloads (Section 4.2). Sec-
by the resulting edge locality, that is, the percentage of uncut
tion 4.3 investigates various parameters of GD.
edges with both endpoints in the same parts. The metric
represents the fraction of local messages in Giraph jobs
4.1 Two-Dimensional Partitioning and corresponds to a possible reduction in communication
Our initial experiments (see Figure 1) and earlier works [18, between the worker machines. Figure 5 reports the results of
29, 33] indicate that two important dimensions for the per- Hash, GD, and BLP on the public dataset. Unsurprisingly, GD
formance of Giraph jobs are the number of vertices and the and BLP outperform the Hash algorithm in the experiment, as
number of edges. For this reason, we specify two weights the latter keeps only k1 of all the edges in the same component.
(1) (2)
for the vertices, wv = 1 and wv = deg(v) for all v ∈ V . The resulting edge locality of GD and BLP are close for the
Recall that our primary goal is to guarantee almost perfect three graphs, though GD typically achieves a higher locality
balance for the two dimensions, as even a single overloaded by 2% − 5%.
 100 Hash BLP GD vertex partitioning edge partitioning vertex+edge partitioning

87.7
40

29.3
75.2

72.3
80 30

67.5
edge locality, %

17.7
61.2

speedup, %
60.5
58.5

14.2
13.3
13.2
20

12.2
52.8

11.1
10.2
60

9.9

8.5
8.1
50

50

50

7.1
42.1

5.6
10

4.9

4.6
-40.1
-17.1

-12.3
-10.6

-53.7

-80.4
-11.5
37.1

-3.1

-3.7
40 0

25.9
23.2
12.5

12.5

12.5
20 -10
-20
0
k=2 k=8 k=2 k=8 k=2 k=8 all rg
e all rg
e all rg
e all rg
e
- sm -la - sm -l a -sm -la -sm -la
LiveJournal Twitter Friendster PR PR CC CC MF MF HP HP

Figure 5: The percentage of local (uncut) edges produced Figure 7: Speedup of Giraph jobs using various partition-
by the three algorithms for the public graphs with k ∈ {2, 8}. ing strategies relative to Hash measured for Page Rank (PR),
Higher values indicate better solutions. GD achieves higher Connected Components (CC), Hypergraph Clustering (HC),
locality in all cases. and Mutual Friends (MF), which are applied on FB-80B
(small ) and FB-400B (large) graphs. Positive values indi-
100 Hash BLP GD
cate improvements, negative ones indicate regressions. Ver-
tex+edge partitioning always results in performance improve-
80 ment.
edge locality, %

52.09
50.07
48.39

60
43.19

42.73
37.16

32.38

29.12
25.42

40
23.03

4.2 Distributed Graph Processing

15.11
13.01

20
6.25

6.25

6.25

In this section we conduct an experimental evaluation of

0.78

0.78

0.78

0 various graph partitioning strategies for speeding up dis-

k=16 k=128 k=16 k=128 k=16 k=128 tributed graph processing. Here we argue and experimen-
FB-3B FB-80B FB-400B tally demonstrate that multi-dimensional balancing is a suit-
able objective for the application. We experiment with four
Figure 6: The percentage of local (uncut) edges produced by graph algorithms implemented in Giraph. Page Rank and
the three algorithms for various subgraphs of the Facebook Connected Components, are popular benchmarks for veri-
friendship graph with k ∈ {16, 128}. Higher values indicate fying the performance of distributed systems. Page Rank
better solutions. GD achieves higher locality in all cases. iteratively propagates vertex ranks through adjacent edges;
our implementation performs 30 iterations for the algorithm.
For the Connected Components algorithm, we use a simple
Figure 6 shows the experiments on the Facebook friendship label propagation technique in which vertices iteratively up-
graphs. Here we use a larger number of parts, k, which more date their labels based on the minimum label of their neigh-
accurately represent the real-world Giraph use case. Again, bors; for our graphs, the process converges after at most 50
Hash produces solutions having the lowest edge localities. In rounds. The other two algorithms, Hypergraph Clustering
fact, over 99% of the edges are cut using the partitioning and Mutual Friends, are applications for large-scale graph
strategy for an instance with a hundred partitions. This is in analytics. The former is used to find a certain clustering of
agreement with our measurements of the typical percentage of the input graph by converting it to a hypergraph. The latter
cross-worker Giraph messages in the production environment. builds a set of features for friend recommendation. Both
On the other hand, we observe a bigger advantage of GD over applications extensively exchange messages between adjacent
BLP; the locality difference is around 10% − 20% for k = 16 vertices, which adds a significant communication overhead.
and 5% − 10% for k = 128. The balanced label propagation Figure 7 depicts the results of our experiment. Since we are
algorithm, BLP, could be configured to produce better results interested in the impact of various partitioning policies on the
by decreasing its cluster size threshold, c. However, this performance of Giraph, we report the relative differences to
results in an imbalanced solution with ε > 0.05 for the the baseline policy, Hash. Here we measure the total runtime
largest instance with k = 128. Hence, we keep the value of of an application using GD as the partitioning strategy in
c = 1024 for all the experiments. The main difference between three modes, vertex partitioning (one-dimensional balance
FB graphs and publicly available graphs is the number of on vertex count), edge partitioning (balance on edge count),
edges. Therefore, we conjecture that the main reason why on and vertex-edge partitioning (two-dimensional balance both
FB graphs GD performs better compared to other algorithms on vertex and edge counts). Every algorithm is applied in
is poor performance of existing local-search based methods two configurations, small and large. The first one uses the
on large graphs in the multi-dimensional case. This is most FB-80B graph and a cluster with 16 worker machines, while
obvious in Figure 6 for k = 128 as one can see that GD is the second one process FB-400B using 128 workers.
gaining a larger advantage over BLP as the size of the graph The key finding is that one-dimensional partitioning cannot
grows (3B → 80B → 400B). However, more experiments provide consistent benefits across all the Giraph applications.
might be required to verify the hypothesis. In fact, we observe performance regression for some instances,
Overall we observe that GD generates solutions of higher in particular, when the number of utilized worker machines
quality than BLP and Hash on the examined instances. There- is large, that is, k = 128. In this scenario, we notice a few
fore, we utilize the algorithm to experiment with distributed workers whose running time is significantly larger than the
graph processing in the next section. average; see Figure 1. Since in Giraph (and other vertex-
Table 2: Impact of partitioning policy on the running time
100 LiveJournal 100 orkut

edge locality, %

edge locality, %
and the amount of sent messages across 128 Giraph workers
for the Page Rank application applied on the FB-400B graph. 90 90
The numbers are average values over 30 iterations. 80 80
70 70
Partitioning Runtime, sec Communication, GB 60 60
mean max stdev mean max stdev 50 0 25 50 75 100 50 0 25 50 75 100
Hash 95 102 27 69.5 69.6 2.4 Iteration Iteration
vertex 93 143 25 18.6 47.6 6.8
edge 82 120 22 25.7 38.2 5.9 Step length: 10 5 2 1
vertex-edge 84 88 21 29.1 30.6 2.8
Figure 8: Comparison of step choices for GD with fixed step
(t) (t+1)
√ that is, kx −x
length, k2 = const, for 100 iterations and
ξ = n/100. Step length 2 · ξ results in good performance.
centric systems) the computation is split into a number of
supersteps that end with a global synchronization barrier,
100 LiveJournal 30 LiveJournal

max imbalance, %
the performance is determined by the slowest worker. Notice

edge locality, %
that a similar phenomena regarding the vertex partitioning 90
has been observed in earlier works [3, 18, 20, 40]. In contrast, 20
80
the two-dimensional partitioning always results in a speedup 10
70
over the default Hash strategy. The improvement is in the
order of 10% − 30% for the examined applications. 60 0
To get a deeper understanding of the source of perfor- 50 0 25 50 75 100 0 25 50 75 100
mance differences, we analyze the detailed logs for the Page Iteration Iteration
Rank application using a cluster with 128 worker machines.
100 orkut 30 orkut

max imbalance, %
edge locality, %
Table 2 shows the measurements of the mean, maximum,
90
and standard deviation of the time to compute a superstep 20
by all the workers. The results indicate that the with the 80
hash-based partitioning, the workers are idling on average for 70 10
7 seconds per superstep waiting for the slowest one to com- 60 0
plete the work. With one-dimensional partitioning the idling 50 0 25 50 75 100 0 25 50 75 100
time is much longer, 50 seconds for vertex-based partitioning Iteration Iteration
and 38 seconds for edge-based one, which is the primary
reason for the performance regression. The two-dimensional nonadaptive adaptive adaptive + vertex fixing
partitioning results in a more even load across the workers
delivering a 13.2% speedup. Table 2 also indicates a signifi- Figure 9: Quality and imbalance comparison of GD 1) with-
cant communication reduction over the baseline partitioning, out adaptive step size, 2) with adaptive step size and 3) with
as measured by the total size of messages sent between the adaptive step size with vertex fixing. The left side shows
workers via network. For the Page Rank application, the av- edge locality and the right side – maximum imbalance over
erage reduction is correlated with the edge locality of the all dimensions. For nonadaptive and adaptive strategies the
corresponding partitioning. However, an unbalanced parti- changes in the number of cut edges and imbalance in the last
tioning causes some workers to use more memory resources iteration are due to fixing in the end of the algorithm the
and become a bottleneck for graph processing. accumulated imbalance resulting from “one-shot” alternating
Finally, we emphasize that the timings analyzed in the projection. Using GD with adaptive step size and vertex fixing
section exclude the running times of the partitioner itself. results in better locality and preserves almost perfect balance
This is realistic for our use case in which the same friendship during algorithm execution.
graph is expected to be processed multiple times for various
tasks. Thus, the extra overhead incurred by a partitioning
strategy is amortized among several runs. between iterations. A natural scaling parameter for the step
√
length is n as it corresponds to the distance between the
4.3 Parameters of GD initial solution x0 = 0 and an integral solution of the form
In this section we perform an experimental comparison of {−1, 1}n . As we show√in Figure 8 for various graphs, a good
n
various choices of the projection step algorithm in GD and choice of step size is 2 100 , where 100 is the limit we set on the
study its convergence properties. Unless specified otherwise, number of iterations due to the constraints on the execution
we use two-dimensional GD in the following setting: 1) bal- running time.
ance is required with respect to the number of vertices and In Figure 9 we show how adaptive step size and vertex
their degrees, 2) in the projection step we use “one-shot” fixing affect the performance of the algorithm. Note that
alternating projection (see Section 3.1), 3) we use adaptive compared with other methods vertex fixing not only im-
step size and vertex fixing as described in Section 3.2. proves quality but also preserves almost perfect balance even
Since behavior of gradient descent algorithms can depend when simple “one-shot” alternating projection is used. Fi-
on selection of the step size parameters, we used experiments nally, in Figure 10 we show analysis of performance of the
to establish convergence of GD with different choices of these algorithm under different choices of the projection method.
parameters. In particular, our implementation aims to ensure The results show that the exact projection algorithm with
that the step length kx(t) −x(t+1) k2 remains close to constant sufficiently large allowed imbalance leads to the best perfor-
 100
edge locality, % LiveJournal 100 orkut the running time of the SHP algorithm exceed the values by

edge locality, %
90 90 a factor of 1.5 − 2 on the same cluster configuration.
80 80
70 70 5. CONCLUSION
60 60 We introduced a new algorithm for the Multi-
50 0 50 0 Dimensional Balanced Graph Partitioning problem,
25 50 75 100 25 50 75 100 which produces balanced partitions according to multiple
Iteration Iteration user-specified weight functions while maintaining high edge
exact ( =0.1) exact ( =0.001) locality. Our computational experiments indicate that this
exact ( =0.01) alternating algorithm is scalable for large real-world graphs, and outper-
forms existing techniques for the problem. Resulting parti-
tions allow one to achieve moderate speedups in computa-
Figure 10: Quality comparison of GD with various projection tional time for several real-world Giraph workloads. This is
methods. We compare exact projection with various allowed in contrast with balancing on just one dimension, that can
imbalance parameters and “one-shot” alternating projection. result in worsen performance.
Allowing more imbalance typically results in higher quality. One of the interesting directions for future work is incorpo-
“One-shot” alternating projection, which we choose as our rating a wider range of balancing requirements, for example,
default implementation option due to its efficiency, produces those that can depend on the resulting partitioning itself,
results comparable with the exact projection. Dykstra’s pro- such as the number of local edges and the maximum number
jection produces the same results as the exact projection, of edges going between any pair of parts in the resulting
and therefore is not shown. partition. The latter quantity can substantially affect perfor-
mance of the distributed computation in Giraph-like systems,
since communication between different machines depends
mance. Larger imbalance permits more partitions, possibly on the number of edges going between the corresponding
including ones with better locality, allowing the overall al- parts. The proposed algorithm cannot directly handle such
gorithm to find results with higher locality. However, the solution-dependent weight functions as they are not express-
alternating projections algorithm can often be used to achieve ible through an a priori fixed collection of weight functions.
similar performance. This is likely due to the fact that the An interesting theoretical question is finding a fast algo-
alternating projections algorithm despite not computing the rithm for exact projection for d > 2. As we show in in the
projection outputs a point close enough to it. full version, it is possible to use nested binary search to find
{ λj } (and therefore the projection) with arbitrary precision.
4.4 Multi-Dimensional Experiments Unfortunately, the running time of the suggested algorithm is
We perform experiments for d = 3 and d = 4 to illustrate unknown, because it is unclear how to estimate left and right
the performance of our algorithms in the multi-dimensional bounds for the binary search. Determining these bounds gives
r −l
case. For these multidimensional experiments in addition an algorithm with running time O(n · di=1 log j δ j ), where
Q
to balancing on the number of vertices and edges we also lj and rj are bounds for λj and δ is the required precision.
balance based on the following additional vertex weights: Another interesting theoretical question is understanding
• Pagerank. We use Pagerank to model activity level of the convergence properties of our algorithm (or a similar
a node. High Pagerank likely means that the vertex is gradient descent based method) under some assumption
accessed often; hence, balancing on the value can be about the spectral properties of the graph. We see this as a
beneficial for load balancing purposes. challenging open problem. While noisy gradient descent is
• Sum of neighbor degrees. We also use the sum of degrees known to have fast convergence for non-convex optimization
over neighbors of a vertex as a weight function. We subject to equality constraints, if inequality constraints are
choose the sum of neighbor degrees as a proxy for the allowed, convergence analysis is open [16].
size of the 2-hop neighborhood of a vertex, which is
computationally expensive to compute for large graphs.
We run experiments on LiveJournal and Orkut graphs,
APPENDIX
as well as on sx-stackoverflow graph (2.6M vertices, 28M A. PROJECTION FOR D = 2
edges) – the largest SNAP graph which is not a social net- In this section we introduce a randomized O(n log n)-time
work. We compare our algorithm with METIS with allowed algorithm for finding projection for d = 2. Recall from Sec-
imbalance of 0.5%. The results are presented in Table 3. tion 2.2 that for y ∈ Rn we need to find λ∗ = (λ∗1 , λ∗2 ) such
They indicate that METIS achieves poor balance for multiple
that h(1) (λ∗ ) = c1 and h(2) (λ∗ ) = c2 . For λ = (λ1 , λ2 ) we
constraints and that GD outperforms METIS by almost all n
(j)
define h(j) (λ) =
P
parameters in most cases. hi (λ) for j ∈ {1, 2}, where
i=1
4.5 Performance Analysis  (j) P (k)
< yi − 1
wi
 if λk wi
Finally, we analyze scalability of our algorithm. Our results 
 k
(j) (k)
 P
hi (λ) = −wi
are obtained on a Hadoop cluster of 128 workers; each of the (j) if λ k wi > yi + 1
machines is a dual-node 2.4 GHz Intel Xeon E5-2680 with  k
 (j) P (k)
wi (yi − λk wi ) otherwise.

256GB RAM. Figure 11 reports the running time of GD in

k
machine-hours on FB-X graphs of various size with balance
on two dimensions. We observe a near-linear growth of the Once we find (λ∗1 , λ∗2 ) we can compute the coordinates of x
(1) (2)
running time with the size of the input graph. In comparison, as xi = [yi − wi λ∗1 − wi λ∗2 ]. We introduce an auxiliary
Table 3: Comparison of GD with METIS. The results show that for high-dimensional balanced partitioning, METIS cannot
guarantee balance across all dimensions. Better results shown in bold.

LiveJournal Orkut sx-stackoverflow

GD METIS GD METIS GD METIS
locality, % 91.71 93.74 88.36 86.52 75.82 80.41
d = 2: balance on vertices imbalance 0.04 0.5 0.02 0.7 0.04 0.6
and degrees memory, MB 2635 4085 4673 10259 1587 4113
time, sec 117 44 203 92 68 55
d = 3: balance on vertices, locality, % 88.74 73.36 89.55 62.1 76.8 60.09
imbalance 0.05 30 0.02 1.6 0.1 6.5
degrees and
memory, MB 2711 4802 4697 12271 1627 4985
sum of neighbor degrees
time, sec 140 66 196 303 76 131
d = 4: balance on vertices, locality, % 87.93 74.36 75.58 65.08 77.04 78.54
degrees, imbalance 0.5 38 2.7 20 0.4 3.8
sum of neighbor degrees memory, MB 2939 4839 4896 12294 1754 5013
and pagerank time, sec 227 66 240 297 88 142

The intuition behind the algorithm is then as follows (in

order to achieve the best performance the exact details differ
machine-hours

120
slightly from this simplified presentation). Suppose we could
80 find a region that contains some solution λ∗ . Then since
constraint functions are linear inside the region, in order to
40
find λ∗ we could solve a system of linear equations over λ1
0 and λ2 . We identify such region, with binary search over
0 200 400 600 800
λ1 by using monotnicity of ∆. We consider only a finite set
number of edges / 10^9 of values: λ1 -coordinates of intersections of boundary lines.
Since there are O(n) boundaries, there are O(n2 ) intersec-
Figure 11: Scalability of the distributed implementation tions(e.g., in Figure 12 we consider only points a, b, c and d).
of GD on FB-X graphs of various size. The results indicate a Hence O(log n) iterations of binary search suffice. The only
near-linear dependence of the running time on the number difference between Algorithm 2 and the above approach is
of edges in the graph. that after the binary search on λ1 we still have to try O(n)
regions to identify the exact region which contains λ∗ (see
Algorithm 2 for the details).
function ∆ which we use to solve the above problem using Now consider one iteration of the binary search. Let λl1 and
r
binary search: 1 λ1 be its current boundaries. Let Λ be a set of all intersection
points (λ1 , λ2 ) such that λ1 ∈ (λl1 , λr1 ). Since ∆ is monotone,
Definition A.1. Suppose that λ1 is such that there exists
for any λ01 we can use binary search by checking whether
λ2 for which the constraint h(2) (λ1 , λ2 ) = c2 is satisfied. Then
λ∗ is greater or less than λ01 through a comparison of ∆(λ01 )
we define ∆(λ1 ) , h(1) (λ1 , λ2 ). and c1 (lines 9-12). Computing ∆(λ01 ) requires solving the
We now describe an O(n log n)-time algorithm for finding one-dimensional problem over λ2 discussed in Section 2.3
(λ∗1 , λ∗2 ). The algorithm is shown as Algorithm 2. It takes as and thus can be done in O(n) time.
a parameter a Boolean variable ∆+ indicating whether ∆ is In order to have binary search run in O(log n) iterations
an increasing or decreasing function. We run the algorithm it suffices to find a value λ01 ∈ (λl1 , λr1 ) which with constant
under both assumptions and select a solution satisfying the probability splits Λ into two subsets of points, those with
constraints. λ1 > λ01 and with λ1 < λ01 respectively, of size at most
2
We outline the main ideas behind Algorithm 2 below. 3
n each. In particular, it suffices to sample a uniformly
Consider the (λ1 , λ2 ) plane partitioned by the following lines random point (λ01 , λ02 ) from Λ. The following lemma bounds
(which we call boundary lines): the overall running time of these sampling steps.
(1) (2)
yi − λ1 wi − λ2 wi =1 Lemma A.1. The overall time required for sampling ran-
(1) (2) dom points from Λ in line 8 of Algorithm 2 is O(n log n).
yi − λ 1 wi − λ 2 wi = −1,
for all i. Let L be the set of boundary lines (line 2). We refer Proof. Consider three cases:
to the subsets of the plane resulting from its partition by the Case 1. |Λ| > n log n. We sample O(n) uniformly random
boundary lines as regions (see Figure 12 where the regions pairs of lines from L and find an intersection of each pair
are referred to as { Ti }). Boundary lines separate the plane (assume no parallel lines which can be handled separately).
into half-planes corresponding to the different cases in the Since the number of lines is O(n) w.h.p. we sample at least
(j) one intersection which lies in Λ. The last condition can be
definitions of the corresponding hi . Therefore, inside each
(j)
checked in O(n) time and if it doesn’t hold then we conclude
region all hi are linear and hence h(j) are also linear. that w.h.p. |Λ| ≤ n log n. We then compute S, the set of all
1
In the full version of the paper we will show that the ∆ is points in Λ in O(n log n) time and proceed to the second
well-defined and monotone. case. We show how to compute Λ in the full version.
 Algorithm 2: Function returning λ∗1 , λ∗2 for given 2- λl1 λ2 λr1
dimensional problem.
(1) (2)
input : { wi }, { wi > 0 }, { yi }, c1 , c2
T7 T8 T9
parameter: ∆+ ∈ {true, f alse} indicating whether ∆ a b
increases
output : (λ∗1 , λ∗2 ) T4 T5 T6
1 Function Project-2D λ1
/* L = set of lines parameterized by
(1) (2) c d
(yi , wi , wi , ±1) corresponding to lines of
(1) (2)
the form yi − λ1 wi − λ2 wi = ±1 */ T1 T2 T3
(1) (2)
2 L := { (yi , wi , wi , ±1) | i ∈ [n] }
3 λl1 := −∞, λr1 := +∞

/* Run binary search */ Figure 12: Example of regions for n = 2, y = (0, 0). Bound-
4 while true do ary lines are λ1 +λ2 = ±1 and λ2 = ±1. These lines partition
5 Λ := set of intersection points (λ1 , λr ) of (λ1 , λ2 )-plane into nine regions T1 , T2 , . . . , T9 . Intersection
lines in L such that λ1 ∈ (λl1 , λr1 ) points are a, b, c and d. Current intersection points consid-
6 if Λ = ∅ then ered by the algorithm (between λl1 and λr1 ) are shown in
7 break red.
8 Sample a uniformly random intersection
point (λ01 , λ02 ) from Λ Theorem A.2. There exists an O(n log n)-time random-
9 if ∆(λ01 ) > c1 then ized algorithm returning λl1 and λr1 such that:
10 If ∆+ set λr1 := λ01 , otherwise set λl1 := λ01 1. No intersections of boundary lines in [λl1 , λr1 ],
11 else 2. There exists a solution (λ†1 , λ†2 ) such that λ†1 ∈ [λl1 , λr1 ].
12 If ∆+ set λl1 := λ01 , otherwise set λr1 := λ01
After we find λl1 and λr1 as in Theorem A.2 we show that
there are only O(n) regions which can contain a solution
13 Let { Rt }Tt=1 be a partition of (λl1 , λr1 ) × R by and we can check them in O(n log n) time. The following
boundary lines, sorted from bottom to top theorem completes the proof of Theorem 1.1 for d = 2:
14 Compute coefficients for the system of
linear equations for R1 (as in Theorem A.3) Theorem A.3. If there exists a solution λ† such that
15 for t = 1 . . . T do λ†1∈ (λl1 ; λr1 ) and no intersection points are between (λl1 ; λr1 )
16 Let (λ1 , λ2 ) be a solution to the linear then λ∗ can be found in O(n log n) time.
system for Rt Proof. We show how to find λ∗ in lines 13-19 of the
17 if (λ1 , λ2 ) ∈ Rt then algorithm. Consider set S = (λl1 , λr1 ) × R. Let { Rt }Tt=1 be
18 return (λ1 , λ2 ) the partition a of S into parts lying between the boundary
lines. Since S doesn’t contain boundary intersections and
19 Update the coefficients corresponding to
there are O(n) boundaries, the size of the partition is O(n).
crossing the boundary line between Rt
For each Rt we solve the following system of equations over
and Rt+1 as shown in Theorem A.3
λ1 and λ2 :
P n
(1)
 hi (λ1 , λ2 ) = c1 ,

i=1
n
Case 2. n ≤ |Λ| ≤ n log n. Note that in this case Λ =  h(2)
P
i (λ1 , λ2 ) = c2

{ (λ1 , λ2 ) ∈ S | λ1 ∈ (λl1 ; λr1 ) }, where S is as defined above. i=1
We sample O(n) random points from S so that w.h.p. we get Since no boundary line crosses Rt , it is a subset of some
at least one point from Λ. As before, if this doesn’t happen, (1) (2)
region. Therefore, hi and hi are linear inside Rt , meaning
we conclude that w.h.p. |Λ| < n and proceed to the last case. that the above system becomes a system of linear equations.
Case 3. |Λ| < n. In this case we maintain Λ directly. If the solution to the system belongs to Rt , then we can take
When we sample a random point (λ01 , λ02 ) ∈ Λ, we remove it as λ∗ . Thus it only remains to show how to find coefficients
from Λ all points on one of the sides from λ01 as directed by for the system in O(n log n) total time.
the binary search. Recall that in Algorithm 2 we assume that { Rt } are
In each of the cases above each of O(log n) iterations takes sorted from bottom to top. For R1 we find the linear system
O(n) time and pre-/post-processing between the cases takes coefficients in O(n) time. Assume that the Rt are already
O(n log n) time. Thus sampling takes O(n log n) time. computed. To find the coefficients for next set Rt+1 , notice
Using the above algorithm we can find λl1 and λr1 such that Rt and Rt+1 are separated by some boundary line. This
(j)
that there are no intersection points between them. Since line corresponds to some hi and therefore crossing it will
(j)
there are O(log n) iterations and each of them requires O(n) change the coefficient of only this hi , and the coefficients
time on average, the total running time is O(n log n). This can be recomputed in O(1) time. Since there are O(n) bound-
completes a proof of the following theorem (corresponding ary lines, the overall time for recomputation is also O(n).
to lines 2-12 of the algorithm). Taking sorting of { Rt } into account, the total running time
is O(n log n).
B. REFERENCES December 5-10, 2016, Barcelona, Spain, pages
2973–2981, 2016.
[1] Dykstra’s projection algorithm. [18] J. E. Gonzalez, Y. Low, H. Gu, D. Bickson, and
[Link] C. Guestrin. Powergraph: distributed graph-parallel
projection_algorithm. Accessed: 2019-02-13. computation on natural graphs. In OSDI, volume 12,
[2] Apache Giraph. [Link] page 2, 2012.
[3] Z. Abbas, V. Kalavri, P. Carbone, and V. Vlassov. [19] J. E. Gonzalez, R. S. Xin, A. Dave, D. Crankshaw,
Streaming graph partitioning: An experimental study. M. J. Franklin, and I. Stoica. GraphX: graph
PVLDB, 11(11):1590–1603, 2018. processing in a distributed dataflow framework. In
[4] A. Amir, J. Ficler, R. Krauthgamer, L. Roditty, and Proceedings of the 11th USENIX conference on
O. S. Shalom. Multiply balanced k -partitioning. In Operating Systems Design and Implementation, pages
LATIN 2014: Theoretical Informatics - 11th Latin 599–613. USENIX Association, 2014.
American Symposium, Montevideo, Uruguay, March 31 [20] Y. Guo, S. Hong, H. Chafi, A. Iosup, and D. Epema.
- April 4, 2014. Proceedings, pages 586–597, 2014. Modeling, analysis, and experimental comparison of
[5] A. Anandkumar and R. Ge. Efficient approaches for streaming graph-partitioning policies. Journal of
escaping higher order saddle points in non-convex Parallel and Distributed Computing, 108:106–121, 2017.
optimization. In Proceedings of the 29th Conference on [21] P. Jain and P. Kar. Non-convex Optimization for
Learning Theory, COLT 2016, New York, USA, June Machine Learning. ArXiv e-prints, Dec. 2017.
23-26, 2016, pages 81–102, 2016. [22] I. Kabiljo, B. Karrer, M. Pundir, S. Pupyrev,
[6] C. Avery. Giraph: Large-scale graph processing A. Shalita, Y. Akhremtsev, and A. Presta. Social hash
infrastructure on Hadoop. Proceedings of the Hadoop partitioner: A scalable distributed hypergraph
Summit. Santa Clara, 11(3):5–9, 2011. partitioner. PVLDB, 10(11):1418–1429, 2017.
[7] K. Aydin, M. Bateni, and V. S. Mirrokni. Distributed [23] G. Karypis and V. Kumar. Metis – unstructured graph
balanced partitioning via linear embedding. In partitioning and sparse matrix ordering system, version
Proceedings of the Ninth ACM International Conference 2.0. Technical report, 1995.
on Web Search and Data Mining, San Francisco, CA, [24] G. Karypis and V. Kumar. Multilevel algorithms for
USA, February 22-25, 2016, pages 387–396, 2016. multi-constraint graph partitioning. In Proceedings of
[8] D. P. Bertsekas. Nonlinear programming. Athena the 1998 ACM/IEEE Conference on Supercomputing,
scientific Belmont, 1999. SC ’98, pages 1–13, Washington, DC, USA, 1998. IEEE
[9] C.-E. Bichot and P. Siarry. Graph partitioning. John Computer Society.
Wiley & Sons, 2013. [25] B. W. Kernighan and S. Lin. An efficient heuristic
[10] S. Boyd and J. Dattorro. Alternating projections. 2003. procedure for partitioning graphs. Bell system technical
[11] S. Boyd and L. Vandenberghe. Convex optimization. journal, 49(2):291–307, 1970.
Cambridge university press, 2004. [26] R. Krauthgamer, J. Naor, and R. Schwartz.
[12] A. Buluç, H. Meyerhenke, I. Safro, P. Sanders, and Partitioning graphs into balanced components. In
C. Schulz. Recent advances in graph partitioning. In Proceedings of the Twentieth Annual ACM-SIAM
Algorithm Engineering - Selected Results and Surveys, Symposium on Discrete Algorithms, SODA 2009, New
pages 117–158. 2016. York, NY, USA, January 4-6, 2009, pages 942–949,
[13] D. Delling, A. V. Goldberg, I. P. Razenshteyn, and 2009.
R. F. F. Werneck. Exact combinatorial [27] J. D. Lee, B. Recht, R. Salakhutdinov, N. Srebro, and
branch-and-bound for graph bisection. In Proceedings J. A. Tropp. Practical large-scale optimization for
of the 14th Meeting on Algorithm Engineering & max-norm regularization. In Advances in Neural
Experiments, ALENEX 2012, The Westin Miyako, Information Processing Systems 23: 24th Annual
Kyoto, Japan, January 16, 2012, pages 30–44, 2012. Conference on Neural Information Processing Systems
[14] L. Dhulipala, I. Kabiljo, B. Karrer, G. Ottaviano, 2010. Proceedings of a meeting held 6-9 December 2010,
S. Pupyrev, and A. Shalita. Compressing graphs and Vancouver, British Columbia, Canada., pages
indexes with recursive graph bisection. In Proceedings 1297–1305, 2010.
of the 22Nd ACM SIGKDD International Conference [28] J. Leskovec and A. Krevl. SNAP Datasets: Stanford
on Knowledge Discovery and Data Mining, KDD ’16, large network dataset collection.
pages 1535–1544, New York, NY, USA, 2016. ACM. [Link] June 2014.
[15] R. L. Dykstra. An algorithm for restricted least squares [29] L. Li, R. Geda, A. B. Hayes, Y. Chen, P. Chaudhari,
regression. Journal of the American Statistical E. Z. Zhang, and M. Szegedy. A simple yet effective
Association, 78(384):837–842, 1983. balanced edge partition model for parallel computing.
[16] R. Ge, F. Huang, C. Jin, and Y. Yuan. Escaping from Proceedings of the ACM on Measurement and Analysis
saddle points - online stochastic gradient for tensor of Computing Systems, 1(1):14, 2017.
decomposition. In Proceedings of The 28th Conference [30] Y. Low, D. Bickson, J. Gonzalez, C. Guestrin,
on Learning Theory, COLT 2015, Paris, France, July A. Kyrola, and J. M. Hellerstein. Distributed graphlab:
3-6, 2015, pages 797–842, 2015. a framework for machine learning and data mining in
[17] R. Ge, J. D. Lee, and T. Ma. Matrix completion has no the cloud. PVLDB, 5(8):716–727, 2012.
spurious local minimum. In Advances in Neural [31] N. Maculan, C. P. Santiago, E. Macambira, and
Information Processing Systems 29: Annual Conference M. Jardim. An O(n) algorithm for projecting a vector
on Neural Information Processing Systems 2016, on the intersection of a hyperplane and a box in r n.
 Journal of optimization theory and applications, [44] S. Wright and J. Nocedal. Numerical optimization.
117(3):553–574, 2003. Springer Science, 35(67-68):7, 1999.
[32] K. Makarychev and Y. Makarychev. Nonuniform graph
partitioning with unrelated weights. In Automata,
Languages, and Programming - 41st International
Colloquium, ICALP 2014, Copenhagen, Denmark, July
8-11, 2014, Proceedings, Part I, pages 812–822, 2014.
[33] C. Martella, D. Logothetis, A. Loukas, and G. Siganos.
Spinner: Scalable graph partitioning in the cloud. In
33rd IEEE International Conference on Data
Engineering, ICDE 2017, San Diego, CA, USA, April
19-22, 2017, pages 1083–1094, 2017.
[34] H. Meyerhenke, P. Sanders, and C. Schulz. Partitioning
complex networks via size-constrained clustering. In
International Symposium on Experimental Algorithms,
pages 351–363. Springer-Verlag New York, Inc., 2014.
[35] J. Nishimura and J. Ugander. Restreaming graph
partitioning: Simple versatile algorithms for advanced
balancing. In Proceedings of the 19th ACM SIGKDD
International Conference on Knowledge Discovery and
Data Mining, pages 1106–1114, New York, NY, USA,
2013. ACM.
[36] S. Ou, K. Yang, and A. Liotta. An adaptive
multi-constraint partitioning algorithm for offloading in
pervasive systems. In 4th IEEE International
Conference on Pervasive Computing and
Communications (PerCom 2006), 13-17 March 2006,
Pisa, Italy, pages 116–125, 2006.
[37] K. Schloegel, G. Karypis, and V. Kumar. A new
algorithm for multi-objective graph partitioning. In
Euro-Par ’99 Parallel Processing, 5th International
Euro-Par Conference, Toulouse, France, August 31 -
September 3, 1999, Proceedings, pages 322–331, 1999.
[38] A. Shalita, B. Karrer, I. Kabiljo, A. Sharma, A. Presta,
A. Adcock, H. Kllapi, and M. Stumm. Social Hash: an
assignment framework for optimizing distributed
systems operations on social networks. In Proceedings
of the 13th Usenix Conference on Networked Systems
Design and Implementation, pages 455–468. USENIX
Association, 2016.
[39] J. Sun, Q. Qu, and J. Wright. When are nonconvex
problems not scary? CoRR, abs/1510.06096, 2015.
[40] J. Sun, H. Vandierendonck, and D. S. Nikolopoulos.
VEBO: A vertex-and edge-balanced ordering heuristic
to load balance parallel graph processing. arXiv
preprint arXiv:1806.06576, 2018.
[41] C. E. Tsourakakis, C. Gkantsidis, B. Radunovic, and
M. Vojnovic. FENNEL: streaming graph partitioning
for massive scale graphs. In Seventh ACM International
Conference on Web Search and Data Mining, WSDM
2014, New York, NY, USA, February 24-28, 2014,
pages 333–342, 2014.
[42] J. Ugander and L. Backstrom. Balanced label
propagation for partitioning massive graphs. In Sixth
ACM International Conference on Web Search and
Data Mining, WSDM 2013, Rome, Italy, February 4-8,
2013, pages 507–516, 2013.
[43] S. Verma, L. M. Leslie, Y. Shin, and I. Gupta. An
experimental comparison of partitioning strategies in
distributed graph processing. PVLDB, 10(5):493–504,
2017.