Practical Parallel Algorithms forGraph Coloring Problems in
Numerical Optimization

Assefaw Hadish Gebremedhin

Thesis submitted in partial fulfillment of therequirements for the degree ofDoctor Scientiarum

Department of Informatics
University of Bergen
NorwayFebruary 2003

Acknowledgements
This work was financially supported by the University of Bergen through a
research fellowship at the Department of Informatics.
I am indebted to many people for their contribution to the successful
completion of this work.
I would first like to thank my advisor Professor Fredrik Manne for his
support and guidance throughout the work with this thesis. Fredrik has
always been encouraging, meticulous in reading my drafts, generous withhis time, and his advices on both academic and non-academic matters havebeen invaluable.
I would also like to thank the other people whom I had the pleasure ofworking with on parts of this thesis. Thanks to Alex Pothen for facilitating
my visit to Old Dominion University, Virginia, USA in the spring of 2001,
for making the stay enjoyable, and for our continued collaboration sincethen. Thanks also to the Pothen family for their hospitality. I thank Jens
Gustedt, Isabelle Guérin Lassous and Jan Arne Telle for our joint research
work and friendship.
Thanks to friends and colleagues at the Department of Informatics forcreating a warm and reliable working place. I also thank all my fellow countrymen living and studying in Bergen for their friendship and the memorable
get-togethers.
I am grateful to my parents Hadish Gebremedhin and Alganesh Tsegay,
my sisters Menen and Wintana, my brothers Theodros, Fasil, Alula andEyob for their support and encouragement.
Finally, I thank my wife Tsigeweini for her love, support and patience.
Without her constant encouragement and never-fading understanding, I
could not have done this.
I dedicate this thesis to our little daughter Melody with love.
Bergen, February 2003
Assefaw Hadish Gebremedhin

1Introduction

Partitioning a set of objects into groups according to certain rules is a fundamental process in mathematics and computer science. Often, the rules that
determine whether or not a pair of objects can belong to the same groupare conceptually simple. Graph coloring is one mathematical abstraction forsuch partitioning rules.
A graph represents binary relationships. It can be visualized as a diagram consisting of points called vertices, and lines called edges. A vertex
represents some physical entity or abstract concept. An edge, which joinsexactly two vertices, represents the relationship between them.
In its standard form, the graph coloring problem deals with assigningcolors to the vertices of a graph such that adjacent vertices receive different
colors and the number of colors used is minimized. This problem and its
variants model a wide range of real-world problems in scientific, engineering and computing applications. For example, timetabling and scheduling,frequency assignment to mobile telephones, register allocation in compilers,printed circuit testing, computation of sparse Jacobian and Hessian matrices, partitioning of tasks in parallel computation, are just a few of the many
areas in which some variant of the graph coloring problem is used as aneffective model.
Among the examples mentioned above, computation of sparse derivative
matrices is an application area considered in this study. The computationof Jacobian and Hessian matrices is a phase required in many numerical
optimization algorithms. Automatic differentiation and finite differencing
are two widely used techniques for computing, or estimating, the entries ofthese matrices. The computation involved in using the techniques can bemade efficient by exploiting matrix structures such as sparsity and symmetry. This leads to several types of partitioning problems that can naturally
be modeled using variants of the graph coloring problem.
This thesis is mainly about developing algorithms for solving different
graph coloring problems. The coloring problems of our concern are all believed to be intractable. Therefore, we focus on finding polynomial-timealgorithms that yield sub-optimal solutions. In particular, we emphasize
greedy algorithms that give usable solutions quickly.
Like many other industrial problems, the graph coloring instances inour cases can be of very large size. One way of handling large-scale problems effectively is by using parallel computers. A parallel computer consistsof several processors that work on different parts of a common problem.
Today, parallel computers consisting of hundreds of processors are commer1

cially available. The basic idea behind parallel computation is to carry out
several tasks simultaneously and thereby reduce execution time. A successful application of this idea requires, among other things, developing parallel algorithms—algorithms that handle the task-division, coordination, andcommunication among the processors. Most of the coloring algorithms developed in this thesis are parallel.
In general, the design and analysis of algorithms, sequential as well as
parallel, presumes the existence of an underlying computation model. Unlike
the sequential case, there is currently no single unifying standard model
for parallel computation. The literature contains an abundant number ofmodels. Among them, the Bulk Synchronous Parallel (BSP) model stands
out for being proposed as a standard bridging model between hardware
and software in parallel computation. The introduction of BSP stimulated
several studies that lead to various adaptations of the model along different
directions. The Coarse Grained Multicomputer (CGM) model is one such
adaptation. Inspired by the BSP and CGM models, in this study, we propose
a new model called the Parallel Resource-Optimal (PRO) model. The PROmodel is proposed in an attempt to further narrow the gap between theory
and practice in developing efficient and scalable parallel algorithms.
In summary, this thesis deals primarily with the design, analysis, andimplementation of practical parallel algorithms for a variety of graph coloring problems. Besides the standard coloring problem, we consider several
specialized variants that arise in the computation of sparse Jacobian andHessian matrices. Secondarily, it is concerned with the development of asimple and effective parallel computation model suitable for theoretical aswell as practical purposes. The thesis consists of the following five papers:
I.

A.H. Gebremedhin and F. Manne. Scalable Parallel Graph Coloring
Algorithms. Concurrency: Practice and Experience. 2000; 12:11311146.

II.

A.H. Gebremedhin, F. Manne and A. Pothen. Graph Coloring in
Optimization Revisited. January 2003. To be submitted.

III.

A.H. Gebremedhin, F. Manne and A. Pothen. Parallel Distance-k
Coloring Algorithms for Numerical Optimization. In B. Monien and
R. Feldmann (Eds.): Proceedings of the 8th International Euro-Par
Conference, Paderborn, Germany, August 2002, volume 2400 of Lecture Notes in Computer Science, pages 912-921. Springer-Verlag 2002.

2

IV.

A.H. Gebremedhin, I.Guérin Lassous, J. Gustedt and J.A. Telle.
Graph Coloring on Coarse Grained Multicomputers. Discrete Applied
Mathematics, to appear.
A preliminary version appeared In U. Brandes and D. Wagner (Eds):
Proceedings of WG 2000, 26th International Workshop on GraphTheoretic Concepts in Computer Science, Konstanz, Germany, June
2000, volume 1928 of Lecture Notes in Computer Science, pages 184195. Springer-Verlag 2000.

V.

A.H. Gebremedhin, I.Guérin Lassous, J. Gustedt and J.A. Telle.
PRO: a Model for Parallel Resource-Optimal Computation. In Proceedings of HPCS’02, 16th Annual International Symposium on High
Performance Computing Systems and Applications, Moncton, Canada,
June 2002, pages 106-113. IEEE Computer Society Press 2002.

The versions of Papers I, IV and V included in this thesis are contentwise identical to the published versions. Paper III is slightly extended fromits published form: it includes more experimental results and proofs thatwere omitted for space considerations. Paper II is planned to be submitted
for publication; unlike the other papers, it is written with a broader audiencein mind.
The topics addressed by Papers I–V can be classified into three groups:• parallel algorithms for the standard graph coloring problem (Paper I),• application of graph coloring in sparse Jacobian and Hessian computation (Papers II and III), and• parallel computation models (Papers IV and V).
The purpose of this common introduction is to highlight the main results in each paper and to put the results in context. It is organized asfollows. In Section 2 we discuss some parallel processing issues relevant forthe subsequent presentation of the papers. In Sections 3–5, we consider thethree topics listed above in their respective order. In particular, in eachsection, a concise summary of the corresponding paper(s) is given. Paper
IV has aspects that make it relevant to both the first and the last topic.
It is, however, presented under the last topic in Section 5. We conclude
the common introduction in Section 6 with some remarks. The body of thethesis, Papers I–V, follows thereafter.

3

2

Multiprocessor parallel computing

The ever increasing need for faster solutions and for solving large scale problems is one of the major driving forces behind parallel computing. The
primary objective here is to achieve increased computational speed by employing a number of processors concurrently.
Concurrent computation can be done at different levels. Our focus is on
using a multiprocessor parallel computer where a number of processors communicate and cooperate to solve a common computational problem. Multiprocessor parallel computation involves three key ingredients: hardware,
software (programming languages, operating systems and compilers), andparallel algorithms [18]. One can use these ingredients as basis for classifying the different kinds of models used in the context of parallel processing.
Specifically, we identify three distinct, and yet closely related, categories ofmodels called parallel architecture, parallel programming, and parallel computation (algorithmic) model. These categories are introduced merely toframe the discussion in the rest of the current section. The discussion focuses on issues relevant to the presentation of the thesis papers.

2.1

Parallel architecture models

Parallel computers are often divided into different architectural classes depending on such issues as control mechanism (SIMD and MIMD computers),
address space organization (shared address space and message-passing architectures) and interconnection network (static and dynamic networks). Moreinformation on each of these classes can be found in books such as [9, 19].

2.2

Parallel programming models

From a programming point of view, shared address space 1 programming(SASP) and message-passing programming (MPP) are perhaps the two mostwidely used models in contemporary parallel computation.
In the SASP model, programs are viewed as a collection of processes accessing a central pool of shared variables. Data exchange among the processors is achieved by reading from, and writing to, the shared variables. This
programming style is naturally suited to the shared-address-space architecture. In the MPP model, programs are viewed as a collection of processes
with private local variables and the ability to exchange data via explicit mes1

The phrases ‘shared address space’ and ‘shared memory’ are often used interchangeably, both in the literature and here in this thesis. There is, however, a difference inmeaning: the former does not require memory to be physically shared.

4

sage passing. In this model, no variables are shared among processors. Eachprocessor uses its local variables for computation, and whenever necessary,it sends data to, or receives data from, other processors. This programming
model fits naturally to distributed-memory computers.
In using the SASP and MPP models, there is usually a trade-off in the
choice to be made. MPP requires considerably more programming effort
than SASP. This is mainly because MPP requires the data structures ina program to be explicitly partitioned whereas SASP does not. Moreover,
SASP supports incremental parallelization of an existing sequential application whereas MPP does not. On the other hand, applications developed
using MPP are likely to be more scalable than SASP based applications.
Currently, OpenMP [23] and MPI [21] are generally considered to be thestandards for writing portable parallel programs using the SASP and MPP
models, respectively. OpenMP, in which the parallel algorithms presented
in Papers I and III are implemented, is a directive-based fork-join modelfor SASP. The fact that OpenMP can be used to write parallel programs
relatively easily makes its usage appealing 2 .
The experimental work reported in Papers I and III is performed on a
Cray Origin 2000 [24], a machine that supports both SASP and MPP. In thismachine, memory is physically distributed among the processors. However,for operational purposes, it behaves as a shared memory computer withthe operating system taking care of maintaining memory consistency. This
enables any processor to use the entire memory of the system.

2.3

Parallel computation models

A computation model is required to serve two major purposes [1]. First, it isused to describe a real entity, namely, a computer. As such, a computation
model attempts to capture the essential features of a machine while ignoring
less important details of its implementation. Second, it is used as a tool for
analyzing problems and expressing algorithms. In this sense, the model is
not necessarily linked to any real computer but rather to an understandingof computation.
In the realm of sequential computation, the Random Access Machine(RAM) is a standard model that has succeeded in achieving both of these
purposes. It serves as an effective model for hardware designers, algorithm
developers and programmers alike. When it comes to parallel computation,
there is no such unifying standard model. This is mainly due to the complex
2

OpenMP dedicated conference series such as EWOMP and WOMPAT could be cited
as indicators for the increasing interest in using OpenMP.

5set of issues inherent in parallel computation.
The natural parallel analogue of RAM, the PRAM, is an overly simplistic model for parallel computation. It is mainly of theoretical interest as it
fails to capture the features of existing parallel computers. The Bulk Synchronous Parallel (BSP) model [26] is a relatively recent model, proposed
to serve as a standard ‘bridging model’ between hardware (machine architecture) and software (algorithm design and programming). As opposed to
the PRAM, parallel algorithms in the BSP model are organized in distinct
computation and communication phases. Moreover, unlike the PRAM, the
BSP model made parallel computation to be coarse grained, a concept wewill make more precise shortly. Later, a variant of the BSP called the Coarse
Grained Multicomputer (CGM) model was proposed [4, 11, 12]. The CGM
model added a strictly coarse grained communication requirement to the
BSP model. The CGM model uses fewer number of parameters than the
BSP model, a feature that simplifies analysis to a certain extent.
In Paper V we propose a new parallel computation model called theParallel Resource-Optimal (PRO) model. The PRO model inherits the advantages offered by the BSP and CGM models. It compromises between
theoretical and practical considerations in the design of optimal and scalable parallel algorithms. Paper V is introduced in Section 5. However, since
the BSP and CGM models are mentioned rather briefly in the paper, we discuss the key attributes of the two models in Sections 2.3.2 and 2.3.3 below.
The two models are presented here as defined in [4].
We proceed in this section by first recapitulating the main features of
the well known PRAM model, which is the model used in Papers I and III
for theoretical runtime analysis. In the same section we include a discussionof a related complexity class. A good introduction to the PRAM model andparallel algorithm design using the model can be found in books such as[1, 16].
2.3.1

PRAM

The PRAM is an idealized parallel computation model. In its standard
form, it consists of an arbitrarily large number of processors and a sharedmemory of unbounded size that is uniformly accessible to all processors.
The processors share a common clock and operate in lockstep, but they may
execute different instructions in each cycle. The PRAM is therefore a model
for a synchronous shared memory MIMD computer. The PRAM can beclassified into different subclasses depending on how simultaneous memory
accesses are handled. For example, a subclass that allows simultaneous read6

access, but not write access, to a single memory location is called concurrent
read exclusive write (CREW) PRAM.
The PRAM is a model for fine-grain computation as it supposes thatthe number of processors can be arbitrarily large. Usually, it is assumed
that the number of processors is a polynomial in the input size. However,practical parallel computation is typically coarse-grain. In particular, on
most existing parallel machines, the number of processors is several ordersof magnitude less than the input size.
Despite its serious deviation from the nature of real parallel computers, the PRAM model is sometimes used for theoretical runtime analysis. A
PRAM-analysis gives an idea on the ‘computational parallelism’ rendered byan algorithm, leaving communication costs aside. Such an analysis should be
supported by empirical results in order to make reasonable conclusions. In
Papers I and III we use the PRAM model for theoretical analyses and support the claims with experimental results. It should, however, be noted thatthe algorithms in Paper I and III are not actually PRAM-algorithms. First,
they do not assume lockstep synchronization. Second, they are developedin a coarse-grain setting.
The class NC On the PRAM model, a parallel algorithm is considered‘efficient’ if its running time is polylogarithmic, i.e., O(log k n) for some fixed
constant k, while the number of processors it uses is polynomial in the
input size n. The class of problems that can be solved by such efficientparallel algorithms is called NC [16]. By simulation, an NC-algorithm canbe converted into a polynomial time sequential algorithm. Thus, the class
NC is included in the class P, i.e., NC ⊆ P. On the other hand, whether or
not P ⊆ NC is an open problem in complexity theory. The general belief
is that P 6⊆ NC, and hence that there are problems in P that do not haveNC-algorithms. The class of P-complete problems consists of the most likely
candidates for such problems. Informally, a problem is said to be P-complete
if an NC-algorithm for it implies that all problems in P have NC-algorithms.
The NC versus P-complete dichotomy on PRAM suffers from several
drawbacks. Most importantly, the dichotomy excludes practical parallelizability in a coarse-grain setting and secondly it does not take work optimality into account. A more detailed discussion of these drawbacks is given inPaper V. Moreover, the experimental results reported in Papers I and III
demonstrate the relevance of developing coarse-grain algorithms.

7

2.3.2The BSP model

A BSP computer is a collection of processor/memory modules connected bya router that can deliver messages in a point to point fashion between theprocessors. A BSP-algorithm is divided into a sequence of supersteps separated by barrier synchronizations. A superstep has distinct computationand communication phases. In a computation phase the processors perform
computation on data that exists locally at the beginning of the superstep.
In a communication phase, data is exchanged among the processors via therouter. The BSP model uses the four parameters, n, p, L, and g. Parametern is the problem size, p is the number of processors, L is the minimum time
between synchronization steps (measured in basic computation units), andg is the ratio of overall system computational capacity (number of computation operations) per unit time divided by the overall system communication
capacity (number of messages of unit size that can be delivered by the router)per unit time.
In a superstep, a processor may send (and receive) at most h messages.
Such a communication pattern is called an h-relation and the basic task of
the router is to realize arbitrary h-relations.
2.3.3

The CGM model

The CGM model is a specialization of the BSP model in which the numberof parameters involved is reduced to just two—p and n—making theoreticalanalysis relatively simpler. The CGM model is a collection of p processors,each with O(n/p) local memory, interconnected by a router that can deliver messages in a point to point fashion. A CGM algorithm consists ofan alternating sequence of computation rounds and communication rounds
separated by barrier synchronizations. A computation round is equivalent
to the computation phase of a superstep in the BSP model. A communication round consists of a single h-relation with h ≤ n/p. In other words, all
the information sent from one processor to another in one communication
round is packed into one long message, thereby minimizing communication
overhead.
The CGM model is used for designing parallel algorithms for coarse grain
systems, in particular, in the case where n  p. Usually, a lower bound on
n/p, e.g. n/p ≥ pδ , for some δ > 0, is required [4]. The value of δ depends
on the nature of the problem for which a parallel algorithm is sought. The
goodness of a CGM parallel algorithm is often measured by the number ofcommunication rounds it uses [4].

8

3Parallel algorithms for standard graph coloring- Paper I

In the simplest case, the graph coloring problem asks for an assignmentof positive integers (called colors) to the vertices of a simple, connected,
undirected graph such that no two adjacent vertices are assigned the samecolor and the total number of colors used is minimized.
The graph coloring problem is NP-hard. The current best known aplog n)2), where n is the number
proximation ratio for the problem is O(n (log(log n)3of vertices in the graph. Moreover, it is known to be not approximablewithin n1/7− for any  > 0 [3].
Despite these rather pessimistic theoretical results, greedy coloring heuristics are often found to be effective in practice. In some applications, theseheuristics find colorings that are within small additive constants of the optimal coloring [6, 17]. The number of colors used by a greedy heuristic depends
on the order in which the vertices are visited and the choice of color made ateach step. If the smallest possible color is chosen at each step, the number ofcolors used by any such sequential greedy heuristic is at most ∆ + 1, where∆ is the maximum degree in the graph. For some applications, this boundcan be quite acceptable, especially if the coloring is obtained quickly.
Such greedy heuristics are inherently sequential and consequently difficult to parallelize. Specifically, the problem of coloring the vertices of agraph in a prespecified order and using the smallest possible color at eachstep is known to be P-complete [13].
Paper I deals with parallelizing greedy heuristics using the shared-addressspace programming model in a coarse-grain setting. Similar previous efforts
focused on using the message-passing programming model on distributedmemory parallel computers. These efforts yielded no speedup [2, 17]. The
algorithms presented in Paper I overcome this shortcoming.
In Paper I we present a new parallel algorithm that colors a graph G =(V, E) using at most ∆ + 1 colors and has an expected parallel runtime
O(|E|/p) on a CREW PRAM for any number of processors p such thatpp ≤ |V |/ 2|E|.
The algorithm consists of three phases. The first two phases are parallel
while the last one is sequential. In the first phase, the vertex set is simply equally partitioned among the p available processors. Each processor
colors its block of vertices paying attention to already colored (local andnon-local) vertices. In this scenario, a pair of adjacent vertices residing on
different processors may be colored concurrently, and possibly get the samecolor, thereby leading to an inconsistency. In the second phase, inconsis9

tencies are detected in parallel. In the final sequential phase, the detectedinconsistencies are rectified.
In the same paper we present a second algorithm that extends this idea
to reduce the number of colors used. The reduction in the number of colors
is obtained by dividing the first phase of the algorithm just described into
two coloring steps. In the improved algorithm, the coloring in the first stepis used to determine an ordering for a coloring in the second step.
Both algorithms are well suited for the shared address space programming model and are implemented using OpenMP. In agreement with thetheoretical analyses, timing results obtained from experiments conductedon the Origin 2000 using a modest number of processors and graphs from finite element methods and eigenvalue computation show that the algorithms
yield speedup.

4

Graph coloring and efficient computation of sparse
derivative matrices

Computing a derivative matrix, i.e. a Jacobian or a Hessian, is a phase required in many numerical optimization algorithms. In large scale problems,
this phase often constitutes an expensive part of the entire computation.
Hence, efficient methods that exploit matrix structure such as sparsity andsymmetry are often required. Finite difference and automatic differentiation (AD) are two widely used methods in the efficient determination of theelements of derivative matrices. In using these methods, the main goal is tominimize the number of function evaluations or AD passes required. The
pursuit of this goal gives rise to several matrix partitioning problems.
Since the early 70’s a considerable amount of research work has been donein the field of sparse derivative matrix computation [5, 6, 7, 8, 10, 14, 15, 20,
22, 25]. Early studies in this field, such as [10, 22, 25], used matrix oriented
approaches. Later, following the pioneering work of Coleman and Moré[6], graph theoretic approaches proved to be advantageous in understanding,analyzing and solving the matrix problems. In particular, graph coloring—in its several variations—was found to be an effective model. However,existing algorithms and software for such coloring problems are developed for
uniprocessor (sequential) computers. Currently, there is a need for parallel
coloring algorithms that can be used to aid the solution of large-scale PDEconstrained optimization problems on parallel computers. Paper II is an
expository work for, and Paper III is a first step towards, addressing thisneed in a flexible manner.

10

Matrix
Nonsym.
Sym.
Nonsym.
Sym.
Nonsym.
Nonsym.
Sym.
Nonsym.
Sym.
Nonsym.

Partition
1d
1d
2d
1d
2d
1d
1d
2d
1d
2d

Schemedirectdirectdirectsubst.
subst.
directdirectdirectsubst.
subst.

Entries
all
all
all
all
all
some
some
some
some
some

Prob.
P1
P2
P3
P4
P5
P6
P7
P8
P9
P10

Formulation
D2 coloring
D 23 coloring
D 32 bicoloring
acyclic coloring
acyclic bicoloring
restricted D2 coloring
restricted D 32 coloring
restricted D 32 bicoloring
not treated
not treated

(II)
[7]
[8]
[5]
[8]
(II)
(II)
(II)

Table 1: Overview of partition/coloring problems considered in Paper II.

4.1

Unifying Framework – Paper II

In the literature, one finds several types of coloring problems characterizing different types of matrix partitioning problems. Understanding thesimilarities and differences among these problems greatly simplifies the development of algorithms—sequential or parallel—for solving them. Paper II
aims at contributing to this end.
The nature of a partitioning problem in the context of efficient computation of sparse derivative matrices depends on the type of the underlying
matrix and the scenario under which the matrix computation is carried out.
Specifically, the particular problem depends on,
1. whether the matrix to be computed is symmetric or nonsymmetric,
2. whether a one-dimensional partition (involving only columns or rows)
or a two-dimensional partition (involving both columns and rows) isused,
3. whether the evaluation scheme employed is direct (solves a diagonal
system) or substitution based (solves a triangular system), and
4. whether all of the nonzero entries of the matrix or only some of themneed to be determined.
The factors outlined above are mutually orthogonal to each other—potentially resulting in 24 different cases. Within our framework, ten ofthese correspond to practically meaningful partitioning problems whereas
the remaining six do not. In Paper II, we use a unified graph theoretic
framework to study eight of these problems; the remaining two are not considered due to their sophistication. The upper part of Table 1 gives an

11overview of the problems considered in the paper and the respective coloring formulations used. The lower part shows the cases not treated in thepaper.
Paper II is a review as well as research paper that integrates known andnew formulations. In the paper we give the necessary background to motivate and formally introduce the matrix partitioning problems, and thendevelop their equivalent graph coloring formulations. Problems P1–P5 havebeen studied previously whereas P6–P8 are investigated for the first time.
The motivation for introducing the new problems and the corresponding restricted coloring formulations is the fact that the computation of a specified
subset of the nonzero entries of a matrix—for instance, for preconditioning
purposes—can be carried out even more efficiently than the computation of
all of the nonzero entries.
In formulating problems P1–P8 as coloring problems, we use a bipartite
graph to represent a nonsymmetric matrix and an adjacency graph to represent a symmetric matrix. In the case of nonsymmetric matrices, bipartite
graph based formulations are shown to be attractive in terms of flexibility and storage space requirement in comparison with known formulations
based on intersection graphs.
The major contributions made in Paper II can be summarized as follows.
• We propose distance-2 coloring as an alternative, more flexible, formulation for problem P1 as compared to the distance-1 coloring formulation of Coleman and Moré [6] for the same problem.
• We expose the interrelationship among the coloring formulations ofproblems P1 through P5 and identify distance-2 coloring as as a genericmodel.
• For problems P3 and P5, based on the previously known relationshipto graph bicoloring [8], we observe a connection to finding a vertex
cover in a graph.
• Using the insight gained from the unified graph theoretic treatment, wedevelop several simple, sequential, heuristic algorithms for the coloring
formulations of problems P1 through P5.
• We formulate the partitioning problems arising in the computation of
a subset of the entries of a matrix (i.e. problems P6–P8) as restrictedcoloring problems.

12

4.2

Parallel Algorithms – Paper III

In Paper III, we develop several shared address space parallel algorithms
for the coloring formulations of problems P1 and P2. The algorithms are
straightforward extensions of the algorithm for the standard (distance-1)
graph coloring problem presented in Paper I. However, in the cases where
the graph to be colored is relatively dense, randomization is used as a supplementary tool to improve scalability. Our PRAM-analyses show that thealgorithms give almost linear speedup for sparse graphs that are large relative to the number of processors. The algorithms are implemented using
OpenMP and results from experiments conducted on the Origin 2000 using various large graphs show that the algorithms indeed yield reasonable
speedup for a modest number of processors.

5

A new parallel computation model

5.1

Graph coloring on the CGM model – Paper IV

The aim of Paper IV is to adapt the shared address space algorithm forthe standard graph coloring problem given in Paper I to the CGM model.
The CGM model used is basically the same as the CGM model described in
Section 2.3.3 except for the fact that the quality of an algorithm is measured
not using communication rounds but rather using an explicit analysis of theoverall computation and communication cost involved. Particularly, a CGM
algorithm is considered efficient if the parallel runtime of the algorithm,
taking both communication and computation costs into account, yields a
speedup (relative to the complexity of the best known sequential algorithm)
that is a linear function in the number of processors used.
In Paper IV we present such an efficient CGM algorithm that colors aconnected graph G = (V, E) using at most ∆ + 1 colors, where ∆ is the
maximum degree in G. The algorithm is given in two variants: randomized
and deterministic. It is shown that on a p-processor CGM model, where|E|/p > p2 , the proposed algorithms have a parallel running time O(|E|/p)
and an overall (computation and communication) cost O(|E|). These bounds
correspond to the average and worst case for the randomized and deterministic versions, respectively.
In comparison with the original algorithm in Paper I, the CGM coloring algorithm in Paper IV is more complicated. Some of the factors that
contribute to the complexity are distributed-memory related issues such asdata distribution, communication overhead, and load balancing. In termsof the underlying algorithmic idea, a key difference between the algorithms
13

of Paper I and IV is that the latter uses recursion.

5.2

The PRO model – Paper V

The notion of ‘efficiency’ of a CGM algorithm as used in Paper IV motivatedus to refine the definition of the model itself. One of the goals was to keep
the number of parameters in the refined model as low as possible and at
the same time make the algorithm design objective more precise. Anothergoal was to find a model that can serve as a design scheme for the algorithm
developer, i.e., a model that identifies the features of a parallel algorithm
that contribute to its practical scalability on different architectures. The
PRO model is partly a result of this effort.
The novel idea behind the PRO model is the focus on relative resource
optimality and the use of granularity as a quality measure. In this model,
the design and analysis of a parallel algorithm is done relative to the complexity of a specific underlying sequential algorithm. The parallel algorithm
is required to be time and space optimal with respect to the complexityof the reference sequential algorithm. As a consequence of its optimality,
a PRO-algorithm yields linear speedup. The quality of a PRO-algorithm
is measured by an attribute of the model called granularity function. Thisfunction, which is expressed in terms of the input size, shows the number ofprocessors that can be employed without sacrificing optimality.
In Paper V, the PRO model is defined formally and it is systematicallycompared with a selection of other models, including BSP and CGM. The
paper also presents PRO-algorithms for matrix multiplication and one-toall broadcast to illustrate how the model is used in algorithm design andanalysis.

6

Conclusion

The first algorithm introduced in Paper I, and later applied in Paper III,
can be seen as a parallelization technique of general interest. The techniqueconsists of (i) breaking up a given problem into p, not necessarily independent, subproblems of almost equal sizes, (ii) solving the p subproblems
concurrently using p processors, and (iii) resolving any inconsistencies that
may result due to the interdependence among the p subproblems sequentially. The technique assumes that resolving an inconsistency can be done
locally. The success of the technique depends on the overall time spent on
detecting and resolving inconsistencies. Papers I and III demonstrate theeffectiveness of the technique in developing shared-address-space parallel al-

14gorithms for various coloring problems. In these papers, we also showed
different ways, including randomization, by which the technique can be enhanced in the context of coloring problems. In Paper IV we extended thebasic idea behind the technique to develop an efficient CGM-algorithm forstandard graph coloring. This algorithm suggests a way for extending theparallelization technique mentioned earlier. Specifically, if the inconsistencies that could result from (ii) can be determined a priori, then (iii) can be
replaced by a recursive application of the technique.
In Paper II we revisited the role of graph coloring in modeling matrix
partitioning problems that arise in numerical estimation of sparse Jacobianand Hessian matrices. We used a unified graph theoretic framework fordealing with the matrix problems in a flexible manner. As a result, wedeveloped several simple and effective sequential heuristic algorithms. In
Paper III, some of these algorithms were parallelized using the shared address space programming model. The work in Papers II and III can serveas basis for further design and implementation of parallel algorithms for thevarious coloring problems.
In Paper V we proposed PRO as a candidate model for the design andanalysis of efficient, scalable and portable parallel algorithms. Verifying its
utility needs further exploration, both theoretical and experimental.
For more specific directions for further work on each of the topics addressed in this thesis, the reader is referred to the concluding remarks given
at the end of each paper.

References[1] S. G. Akl. Parallel Computation. Prentice Hall, New Jersey, USA, 1997.
[2] J. R. Allwright, R. Bordawekar, P. D. Coddington, K. Dincer, and C. L.
Martin. A comparison of parallel graph coloring algorithms. Technical
Report Tech. Rep. SCCS-666, Northeast Parallel Architecture Center,
Syracuse University, 1995.
[3] G. Ausiello, P. Crescenzi, G. Gambosi, V. Kann, A. MarchettiSpaccamela, and M. Protasi. Complexity and Approximation. Springer,
Berlin, Germany, 1999.
[4] E. Caceres, F. Dehne, A. Ferreira, P. Locchini, I. Rieping, A. Roncato,
N. Santoro, and S. W. Song. Efficient parallel graph algorithms forcoarse grained multicomputers and BSP. In The 24th International

15

Colloquium on Automata Languages and Programming, volume 1256 of
LNCS, pages 390–400. Springer Verlag, 1997.
[5] T. F. Coleman and J. Cai. The cyclic coloring problem and estimationof sparse Hessian matrices. SIAM J. Alg. Disc. Meth., 7(2):221–235,
April 1986.
[6] T. F. Coleman and J. J. Moré. Estimation of sparse Jacobian matrices
and graph coloring problems. SIAM J. Numer. Anal., 20(1):187–209,
February 1983.
[7] T. F. Coleman and J. J. Moré. Estimation of sparse Hessian matrices
and graph coloring problems. Math. Program., 28:243–270, 1984.
[8] T. F. Coleman and A. Verma. The efficient computation of sparse Jacobian matrices using automatic differentiation. SIAM J. Sci. Comput.,
19(4):1210–1233, July 1998.
[9] D. E. Culler, J. P. Singh, and A. Gupta. Parallel Computer Architecture. Morgan Kaufmann Publishers, San Francisco, California, 1999.
[10] A. R. Curtis, M. J. D. Powell, and J. K. Reid. On the estimation ofsparse Jacobian matrices. J. Inst. Math. Appl., 13:117–119, 1974.
[11] F. Dehne. Coarse grained parallel algorithms. Algorithmica Special
Issue on “Coarse grained parallel algorithms”, 24(3/4):173–176, 1999.
[12] F. Dehne, A. Fabri, and A. Rau-Chaplin. Scalable parallel computational geometry for coarse grained multicomputers. International Journal on Computational Geometry, 6(3):379–400, 1996.
[13] R. Greenlaw, H.J. Hoover, and W. L. Ruzzo. Limits to Parallel Computation: P-Completeness Theory. Oxford University Press, New York,USA, 1995.
[14] S. Hossain and T. Steihaug. Computing a sparse Jacobian matrix byrows and columns. Optimization Methods and Software, 10:33–48, 1998.
[15] S. Hossain and T. Steihaug. Reducing the number of AD passes for computing a sparse Jacobian matrix. In G. Corliss, C. Faure, A. Griewank,
L. Hascoët, and U. Naumann, editors, Authomatic Differentiation of Algorithms: From Simulation to Optimization, pages 263 – 270. Springer,
2002.
[16] J. Jájá. An Introduction to Parallel Algorithms. Addison-Wesley, 1992.
16

[17] M. T. Jones and P. E. Plassmann. A parallel graph coloring heuristic.
SIAM J. Sci. Comput., 14(3):654–669, May 1993.
[18] H. F. Jordan and G. Alaghband. Fundamentals of Parallel Processing.
Prentice Hall, New Jersey, USA, 2003.
[19] V. Kumar, A. Grama, A. Gupta, and G. Karypis. Introduction to
Parallel Computing. The Benjamin/Cummings Publishing Company,
Inc., California, 1994.
[20] S. T. McCormick. Optimal approximation of sparse Hessians and its
equivalence to a graph coloring problem. Math. Program., 26:153–171,
1983.
[21] MPI. Message-passing interface standard. http://www.mpi-forum.org/.
[22] G. N. Newsam and J. D. Ramsdell. Estimation of sparse Jacobianmatrices. SIAM J. Alg. Disc. Meth., 4:404–418, 1983.
[23] OpenMP. A proposed industry standard API for shared memory programming. http://www.openmp.org/.
[24] Parallab.
High
Performance
http://www.parallab.uib.no/.

ComputingLaboratory.

[25] M. J. D. Powell and PH. L. Toint. On the estimation of sparse Hessian
matrices. SIAM J. Numer. Anal., 16(6):1060–1074, December 1979.
[26] L. G. Valiant. A bridging model for parallel computation. Communications of the ACM, 33(8):103–111, 1990.

17

Scalable parallel graph coloring algorithms
Assefaw Hadish Gebremedhin∗

Fredrik Manne

Abstract
Finding a good graph coloring quickly is often a crucial phase in
the development of efficient, parallel algorithms for many scientificand engineering applications. In this paper we consider the problem
of solving the graph coloring problem itself in parallel. We present a
simple and fast parallel graph coloring heuristic that is well suited forshared memory programming and yields an almost linear speedup on
the PRAM model. We also present a second heuristic that improveson the number of colors used. The heuristics have been implemented
using OpenMP. Experiments conducted on an SGI Cray Origin 2000
super computer using very large graphs from finite element methods
and eigenvalue computations validate the theoretical run-time analysis.
Key words: graph coloring; parallel algorithms; shared memory
programming; OpenMP∗

Department of Informatics, University of Bergen, N-5020 Bergen, Norway. {assefaw,
Fredrik.Manne}@ii.uib.no1

1

Introduction

The graph coloring problem (GCP) deals with assigning labels (called colors)
to the vertices of a graph such that adjacent vertices do not get the samecolor. The primary objective is to minimize the number of colors used. The
GCP arises in a number of scientific computing and engineering applications.
Examples include, among others, time tabling and scheduling [14], frequency
assignment [6], register allocation [3], printed circuit testing [8], parallel
numerical computation [1], and optimization [4]. Coloring a general graphwith the minimum number of colors is known to be an NP-hard problem [7],thus one often relies on heuristics to obtain a usable solution.
In a parallel application a graph coloring is usually performed in orderto partition the work associated with the vertices into independent subtasks
such that the subtasks can be performed concurrently. Depending on theamount of work associated with each vertex, there are basically two coloring
strategies one can pursuit. In the first strategy the emphasis is on minimizing the number of colors whereas in the second the focus is on speed.
Ascertaining which is more appropriate depends on the underlying problem
one is trying to solve.
If the task associated with each vertex is computationally expensive
then it is crucial to use as few colors as possible. There exist several timeconsuming (iterative) local improvement heuristics for addressing this need.
Some of these heuristics have been shown to be parallelizable [14].
If, on the other hand, the task associated with each vertex is fairly small
and one repeatedly has to find new graph colorings then the overall timeto perform the colorings might take up a significant portion of the entire
computation. See [16] for an example of this case. In such a setting it ismore important to compute a usable coloring fast than spending time on
reducing the number of colors. For this purpose there exist several linear
time, or close to linear time, sequential greedy coloring heuristics. These
heuristics have been found to be effective in coloring graphs that arise froma number of applications [4, 11]. Because of their inherent sequential nature,
however, these heuristics are difficult to parallelize.
This paper focuses mainly on the latter strategy where the goal is todevelop scalable parallel coloring heuristics based on greedy methods. Previous work on developing such algorithms has been performed on distributed
memory computers using explicit message-passing. The speedup obtained
from these efforts has been discouraging [1]. The main justification for using
these algorithms has been access to more memory and thus the ability tosolve problems with very large graphs. It is to be noted that the current
2

availability of shared memory computers where the entire memory can beaccessed by any processor makes this argument less significant now.
The development of shared memory computers has been accompaniedby the emergence of new shared memory programming paradigms of which
OpenMP has become one of the most successful and widely used [17]. OpenMP
is a directive-based, fork-join model for shared memory parallelism.
In this paper we present a fast and scalable parallel graph coloring algorithm suitable for the shared memory programming model. In our context,
scalability of a parallel algorithm is a measure of its capacity to increase
speedup as the number of processors is increased for a given problem size.
Our algorithm is based on first performing a parallel pseudo-coloring of thegraph. The prefix “pseudo” is used to reflect that the coloring might containadjacent vertices that are colored with the same color. To remedy this we
perform a second parallel step where any inconsistencies in the coloring aredetected. These are then resolved in a final sequential step. An analysis
on the PRAM model using p processors for a graph with n vertices and medges shows that the expected number of conflicts from the first stage is
low and for p ≤ √n2m the algorithm is expected to provide a nearly linear
speedup. We also extend this idea and present a second parallel algorithm
that potentially uses fewer colors.
The presented algorithms have been implemented in Fortran90 using
OpenMP on a Cray Origin 2000 super computer. Experimental results on
a number of very large graphs show that the algorithms yield good speedup
and produce colorings of comparable quality to that of their sequential counterparts. The fact that we are using OpenMP makes our implementation
significantly simpler and easier to verify than if we had used a distributed
memory programming environment such as MPI.
The rest of this paper is organized as follows. In Section 2 we givesome background on the graph coloring problem and previous efforts madeto solve it in parallel. In Section 3 we describe our new parallel graphcoloring algorithms and analyze their performance on the PRAM model.
Synchronization overhead and OpenMP issues related to our implementation
are discussed in Section 4. In Section 5 we present and discuss results fromexperiments performed on the Cray Origin 2000. Finally, in Section 6 wegive concluding remarks.

2

BackgroundIn this section we give a brief overview of previous work done on the development of fast coloring heuristics, both sequential and parallel. We begin3

by introducing some graph notation used in this paper.
For a graph G = (V, E), we denote |V | by n, |E| by m, and the degreeof a vertex vi by deg(vi ). Moreover, the maximum, minimum, and average
degree in a graph are denoted by ∆, δ, and δ respectively.
As mentioned in Section 1 there exist several fast sequential coloringheuristics that are very effective in practice. These algorithms are all basedon the same general greedy framework: a vertex is selected according tosome predefined criterion and then colored with the smallest valid color.
The selection and coloring continues until all the vertices in the graph arecolored.
Some of the suggested coloring heuristics under this general framework
include Largest-Degree-First-Ordering (LFO) [18], Incidence-Degree-Ordering(IDO) [4], and Saturation-Degree-Ordering (SDO) [2]. These heuristics
choose at each step a vertex v from the set of uncolored vertices with the
maximum “degree”. In LFO, the standard definition of degree of a vertex
is used. In IDO, incidence degree is defined as the number of already colored adjacent vertices, whereas in SDO one only considers the number ofdifferently colored adjacent vertices. First Fit (FF) is yet another, simple
variant of the general greedy framework. In FF, the next vertex from somearbitrary ordering is chosen and colored. In terms of quality of coloring,these heuristics can in most cases be ranked in an increasing order as FF,
LFO, IDO, and SDO. Note that for a graph G the number of colors used byany sequential greedy algorithm is bounded from above by ∆ + 1. On the
average, however, it has been shown that for random graphs FF is expected
to use no more than 2χ(G) colors, where χ(G) is the chromatic number 1of G [10]. In terms of run time, FF is clearly O(m), LFO and IDO can be
implemented to run in O(m), and SDO in O(n 2 ) [2, 11].
When it comes to parallel graph coloring, a number of the existing fastheuristics are based on the observation that an independent set of vertices
can be colored in parallel. A general parallel coloring scheme based on thisobservation is outlined in Figure 1.
Depending on how the independent set is chosen and colored, Scheme 1
specializes into a number of variants. The Parallel Maximal Independent set(PMIS) coloring is one variant. This is a heuristic based on Luby’s maximal
independent set finding algorithm [15]. Other variants are the asynchronous
parallel heuristic by Jones and Plassmann (JP) [11], and the Largest-DegreeFirst(LDF) heuristic developed independently by Gjertsen Jr. et al. [12] andAllwright et al. [1].
1

The chromatic number of a graph is the optimal number of colors required to color it.

4

Scheme 1
P arallelColoring(G = (V, E))
begin
U ←V
G0 ← G
while (G0 is not empty) do in parallel
Find an independent set I in G0
Color the vertices in I
U ←U \IG0 ← graph induced by U
end-while
end
Figure 1: A parallel coloring heuristic
Allwright et al. made an experimental, comparative study by implementing the PMIS, JP, and LDF coloring algorithms on both SIMD and MIMD
parallel architectures [1]. They report that they did not get speedup for anyof these algorithms.
Jones and Plassmann [11] do not report on obtaining speedup for theiralgorithms either. They state that “the running time of the heuristic is onlya slowly increasing function of the number of processors used”.

3

Block Partition Based Coloring Heuristics

In this section we present two new parallel graph coloring heuristics and give
their performance analysis on the PRAM model. Our heuristics are basedon dividing the vertex set of the graph into p successive blocks of equal size.
We call this a block partitioning. We assume that the vertices are listed in arandom order and thus no effort is made to minimize the number of crossing
edges. A crossing edge is an edge whose end points end up in two different
blocks. Obviously, because of the existence of crossing edges, the coloring
subproblems defined by each block are not independent.

3.1

The First Algorithm

The strategy we employ consists of three phases. In the first phase, theinput vertex set V of the graph G = (V, E) is partitioned into p blocks
as {V1 , V2 , . . . , Vp } such that |Vi | = n/p, 1 ≤ i ≤ p. The vertices in eachblock are then colored in parallel using p processors. The parallel coloring5

comprises of n/p parallel steps with synchronization barriers at the end ofeach step. When coloring a vertex, all its previously colored neighbors, boththe local ones and those found on other blocks, are taken into account. In
doing so, two processors may simultaneously attempt to color vertices thatare adjacent to each other. If these vertices are given the same color, theresulting coloring becomes invalid and hence we call the coloring obtained
a pseudo coloring. In the second phase, each processor p i checks whether
vertices in Vi are assigned valid colors by comparing the color of a vertex
against all its neighbors that were colored at the same parallel step in the firstphase. This checking step is also done in parallel. If a conflict is discovered,
one of the end points of the edge in conflict is stored in a table. Finally,in the third phase, the vertices stored in this table are colored sequentially.
Algorithm 1 provides the details of this strategy and is given in Figure 2.
Algorithm 1
BlockP artitionBasedColoring(G, p)
begin
1. Partition V into p equal blocks V1 . . . Vp , where b np c ≤ |Vi | ≤ d np e
for i = 1 to p do in parallelfor each vj ∈ Vi do
assign the smallest legal color to vertex v jbarrier synchronize
end-forend-for
2. for i = 1 to p do in parallelfor each vj ∈ Vi dofor each neighbor u of vj that has been colored at the same
parallel step do
if color(vj ) = color(u) thenstore min {u, vj } in table A
end-ifend-forend-forend-for3. Color the vertices in A sequentiallyend

Figure 2: Block partition based coloring63.1.1Analysis

Our analysis is based on the PRAM model. Without loss of generality weassume that n/p, the number of vertices per processor, is an integer. Let
the vertices on each processor be numbered from 1 to n/p and the paralleltime used for coloring be divided into n/p time slots. The processors aresynchronized at the end of each time unit t j . This means, at each time unit
tj , processor pi colors vertex vj ∈ Vi , 1 ≤ j ≤ n/p and 1 ≤ i ≤ p.
Our first result gives an upper bound on the expected number of conflicts(denoted by K) created at the end of Phase 1 of Algorithm 1 for a graph inwhich the vertices are listed in a randomly permuted order.
Lemma 3.1 The expected number of conflicts created at the end of Phasen( n−1) ≈ δ(p−1).
1 of Algorithm 1 is at most δ(p−1)
22Proof : Consider a vertex x ∈ V that is colored at time unit t j , 1 ≤j ≤ n/p. Since the neighbors of x are randomly distributed, the expected
number of neighbors of x that are concurrently colored at time unit t j isgiven byp−1
deg(x)(1)
n−1
If we sum (1) over all vertices in G we count each pair of adjacent vertices
that are colored simultaneously twice. Moreover each term in the sum represents only a potential conflict since two adjacent vertices could be coloredsimultaneously and yet be assigned different colors. The sum thus gives anupper bound on the expected number of conflicts. Therefore, we have
E[K] ≤ (1/2)

X p−1
x∈V

n−1

deg(x)

p−1
(2m)
n−1
= (1/2)δ(p − 1)(n/n − 1)
= (1/2)

P

(2)
(3)
(4)

deg(v)In going from (3) to (4), the identity δ = v∈Vn
= 2mn is used. Notethat for large values of n, (n/n − 1) ≈ 1 and the result follows.
2
We now look at the expected run time2 of Algorithm 1. To do so, we
2
introduce a graph attribute called relative sparsity r, defined as r = nm .
Note that 1/r, the ratio of the actual number of edges to the total possible
number of edges, shows the density of the graph. The following lemma states
2

We use the prefix E to identify expected time complexity expressions.

7qthat for bounded degree graphs and for p ≤ 2r , Algorithm 1 provides analmost linear speedup compared to the sequential First Fit algorithm.
Lemma 3.2 On a CREW PRAM, Algorithm
q 1 colors the input graph conr
sistently in EO(∆n/p) time when p ≤2 and in EO(∆δp) time when
p>

q

r
2.

Proof : Note first that since Phase 3 resolves all the conflicts that are
discovered in Phase 2, the coloring at the end of Phase 3 is a valid one.
Both Phase 1 and 2 require concurrent read capability and thus the required
PRAM is CREW. The overall time required by Algorithm 1 is T = T 1 +T2 +T3 , where Ti is the time required by Phase i. Both Phase 1 and 2 consistof n/p parallel steps. The number of operations in each parallel step is
proportional to the degree of the vertex under investigation which is boundedfrom above by ∆. Thus, T1 = T2 = O(∆n/p). The time required by thesequential step (Phase 3) is T3 = O(∆K) where K is the number of conflicts
discovered in Phase 2. From Lemma 3.1, E[K] = O(δp). Substituting yields
T = T1 + T2 + T3 = EO(∆(n/p + δp))(5)

The overall time T is thus determined by how
q n/p compares with δp.

Using the identity δ = 2m
n , we see that for p ≤

n2
2m =

dominates giving an overall running time of EO(∆n/p).

qr
2 , the term n/p
q
For p > 2r , theterm δp dominates and the overall time becomes EO(∆δp).
2
For most practical applications and currently available parallel computers we expect that both δ  n and p  n implying that n > p 2 δ and thus
giving an overall time complexity of O(∆n/p) for Algorithm 1.
The number of colors used by Algorithm 1 is bounded from above by∆ + 1. This follows since the conflicts that arise in Phase 1 are resolved
sequentially. However, we note that there exist instances where the coloring produced by Algorithm 1 can be arbitrarily worse than that of thesequential FF algorithm. To see this, consider a complete bipartite graph
G = (V1 , V2 , E) with |V1 | = |V2 | = n/2, where the vertices in V1 are ordered
before the vertices in V2 and with p = 2. For this setting Algorithm 1 will
use n2 + 1 colors while sequential FF will color the graph optimally using 2
colors.

8

3.2

The Second Algorithm

In this section we show how Algorithm 1 can be modified to use fewer colors.
Our method is motivated by the idea behind Culberson’s Iterated Greedy
coloring heuristic (IG) [5]. IG is based on the following result, stated here
without proof.
Lemma 3.3 (Culberson) Let C be a k-coloring of a graph G, and πa permutation of the vertices such that if C(v π(i) ) = C(vπ(l) ) = c, thenC(vπ(j) ) = c for i < j < l. Then, applying the First Fit algorithm to G
where the vertices have been ordered by π will produce a coloring using k or
fewer colors.
From Lemma 3.3, we see that if FF is re-applied on a graph where thevertex set is ordered such that vertices belonging to the same color class 3 inthe previous coloring are listed consecutively, the new coloring is better orat least as good as the previous coloring. There are many ways in which thevertices of a graph can be arranged satisfying the condition of Lemma 3.3.
One such ordering is the reverse color class ordering [5]. In this ordering,the color classes are listed in reverse order of their introduction. This has
a potential for reducing the number of colors used since one now proceedsby first coloring vertices that could not be colored with low values in theprevious coloring.
Our improved coloring heuristic uses Lemma 3.3 and consists of 4 phases,
one more phase than Algorithm 1. The first phase is the same as Phase 1of Algorithm 1. Let the coloring number used by this phase be ColN um.
During the second phase, the pseudo coloring of the first phase is used to geta reverse color class ordering of the vertices. The second phase consists of
ColN um steps. In each step k, the vertices of color class ColN um−k−1 arecolored afresh in parallel in a similar manner as in Phase 1. The remaining
two phases are the same as Phases 2 and 3 of Algorithm 1. The method justdescribed (Algorithm 2) is outlined in Figure 3.
Each color class at the end of Phase 1 is a pseudo independent set. In
particular any edge within a color class results from a “conflict” edge from
Phase 1. Hence a new block partitioning of the vertices of each color class
results in only a few crossing edges. In other words, the number of conflicts
expected at the end of Phase 2 (K2 ) should be much smaller than the numberof conflicts at the end of Phase 1 (K1 ). Thus, in addition to improving thequality of the coloring, Phase 2 should also provide a significant reductionin the number of conflicts. Note that conflict checking and removal steps
3

Vertices of the same color constitute a color class.

9

Algorithm 2
ImprovedBlockP artitionBasedColoring(G, p)
begin
1. As Phase 1 of Algorithm 1{At this point we have the pseudo independent
sets ColorClass(1) . . . ColorClass(ColNum) }2. for k = ColN um down to 1 do
Partition ColorClass(k) into p equal blocks V 10 . . . Vp0
for i = 1 to p do in parallelfor each vj ∈ Vi0 do
assign the smallest legal color to vertex v j
end-forend-forend-for
3. As Phase 2 of Algorithm 1
4. As Phase 3 of Algorithm 1
end

Figure 3: Modified block partition based coloringare included in Phases 3 and 4 to ensure that any remaining conflicts areresolved.
The following result gives a bound on the expected number of conflicts
at the end of Phase 2 of Algorithm 2.
Lemma 3.4 The expected number of conflicts created at the end of Phase2
2 2≈ 2p n∆ .
2 of Algorithm 2 is at most 2p ∆(∆+1)
n
Proof : From Lemma 3.1, the expected number of conflicts at the end of
Phase 1 is approximately bounded by δp/2. Noting that there are m edges
in the input graph G to Algorithm 2, at the end of Phase 1, the probability
that an arbitrary edge in G is in “conflict” is expected to be no more thanδp2m = p/n. Now consider a color class w from the coloring obtained at theend of Phase 1 of Algorithm 2. Let G0 = (V 0 , E 0 ) be the graph induced by
the vertices of this color class and let n 0 = |V 0 |, m0 = |E 0 |. Further, let x bea vertex in G0 and deg 0 (x) its degree. From the above discussion, we expectthat deg 0 (x) ≤ np deg(x). Using the same argument as in Lemma 3.1, theexpected number of neighbors of x that are concurrently colored at time unit
0tj , for 1 ≤ j ≤ n0 /p, is np−1
0 −1 deg (x). Thus the number of conflicts created10due to the vertices of color class w (denoted by K 0 ) is bounded as follows.
E[K 0 ] ≤X p(p − 1)
x∈V 0

=
≤

deg(x)(6)x∈V 0 deg(x)n0 − 1(7)n(n0 − 1)p(p − 1)
n
2p2 ∆
n

P

(8)

Recall that there are at most ∆ + 1 colors at the end of Phase 1. Therefore,
K2 , the total number of expected conflicts at the end of Phase 2 , isE[K2 ] ≤2p2 ∆(∆ + 1)
n

(9)
2

√
2
2p
) .
Noting that ∆(∆ + 1) ≈ ∆2 , (9) can be rewritten as E[K2 ] ≤ ( √n/∆√√This indicates that if 2p < n/∆, the expected number of conflicts at theend of Phase 2 is less than 1.

4

Implementation Issues

In this section we address the problem of synchronization overhead andillustrate how OpenMP is used in our Fortran90 implementations.

4.1

Synchronization Overhead

The barrier synchronization in Phase 1 of Algorithm 1 is introduced to
identify the parallel step tj (1 ≤ j ≤ n/p) during which a vertex is colored. This information is used for two purposes: (i) in Phase 1 to identify
already colored neighbors of a vertex, and (ii) in Phase 2 to identify theneighbors of a vertex that are colored at the same parallel step as itself.
Although the barrier enables us to realize these purposes, its implementation typically incurs an undesirable large overhead. To overcome this wehave developed an asynchronous version of Algorithm 1 (and consequentlyof Algorithm 2). In the asynchronous version we consider all the neighborsof a vertex under investigation, irrespective of the parallel step during which
they are colored. This is done first when determining the color of a vertex
and then when checking for consistency of coloring. We have implemented
and tested both the asynchronous and synchronous versions of Algorithm
1. In the synchronous version, an OpenMP library routine was utilized torealize barrier synchronization. The obtained results show that the asynchronous version runs faster by a factor of 3 to 5. The relative slow-down
11

!$omp parallel do schedule(static, Bsize)
private(i) shared(vertex)
do i = 1, number of vertices
call assign color synch(vertex(i))
call mp barrier
enddo!$omp parallel do schedule(static, Bsize)
private(i) shared(vertex)
do i = 1, number of vertices
call assign color asynch(vertex(i))
enddo

Figure 4: OpenMP-sketch of Phase 1 of Algorithm 1, synchronous (left) andasynchronous (right)
factor of the synchronous version depends on how often the OpenMP barrier
routine is called. Particularly for a given graph, the relative time spent on
synchronization increases with the number of processors.

4.2

OpenMP

Figure 4 provides a sketch of Phase 1 of both the synchronous and asynchronous versions of Algorithm 1. Here the vertices are stored in the integer
array vertex and the number of vertices per processor is stored in the variable
Bsize. The routines assign color synch(i) and assign color asynch(i), synchronous and asynchronous versions respectively, assign the smallest valid
color to vertex i. It should be noted that even though the synchronousalgorithm visits fewer vertices both when determining the color of a vertex and when checking for consistency, it incurs an extra initial overhead
since one must determine for each of value of p which vertices to check.
The routine mp barrier is an OpenMP library routine that enables barrier
synchronization.
In addition to the standard OpenMP directives we have used data distribution directives provided by SGI to ensure that most cache misses are
satisfied from local memory.
We disallow access to a memory location while it is being written byusing the ATOMIC directive in OpenMP. This makes accessing the color ofany vertex without reading garbage values possible in Phase 1.

5

Experimental Results

In this section, we experimentally demonstrate the performance of the asynchronous versions of Algorithms 1 and 2. In Section 5.1 we introduce the test
graphs used in the experiments and in Section 5.2 we present and discuss
the experimental results. The experiments have been performed on a Cray
Origin 2000, a CC-NUMA machine consisting of 128 MIPS R10000 processors. The algorithms have been implemented in Fortran90 and parallelized

12

Set
Set I
Set I
Set II
Set II
Set III
Set III

Problem
mrng2
mrng3
598a
m14b
dense1dense2n
1,017,253
4,039,160
110,971
214,765
19,703
218,849m
2,015,714
8,016,848
741,934
1,679,018
3,048,477
121,118,458

∆
4
4
26
40
504
1,640

δ
2
2
5
4
116
332

δ
4
4
13
16
309
1,107

pr

2

506
1008
91
117
8
14

χF F
5
5
11
13
122
377

χIDO
5
5
9
10
122
376

Table 1: Test Graphs.
using OpenMP[17]. We have also implemented the sequential versions of FF
and IDO to use as benchmarks.
In these experiments, the block partitioning is based on the ordering of
the vertices as provided in the input graph. In other words, no random
permutation is done on the ordering of the vertices prior to partitioning.

5.1

Test Graphs

The test graphs used in our experiments are divided into three categories
as Problem Set I, II, and III (see Table 1). Problem Sets I and II consistof graphs (matrices) that arise from finite element methods [13]. Problem
Set III consists of matrices that arise in eigenvalue computations [16]. In
addition to providing some statistics about the structure of the test graphs,
Table 1 also lists the number of colors required when coloring the graphs
using our sequential FF and IDO implementations (shown under columns
χF F and χIDO , respectively).

5.2

Discussion

Algorithm 1. Table 2 lists results obtained using the asynchronous version of Algorithm 1. The number of blocks (processors) is given in column
p. Columns χ1 and χ3 give the number of colors used at the end of Phases
1 and 3, respectively. The number of conflicts that arise in Phase 1 arelisted under the column labeled K. The column labeled δ(p−1)
gives the2
theoretically expected upper bound on the number of conflicts as predictedby Lemma 3.1. The time in milliseconds required by the different phases arelisted under T1 , T2 , T3 , and the last column Ttot gives the total time used.
The column labeled Spar lists the speedup obtained compared to the timetot (1)
used by running Algorithm 1 on one processor (S par (p) = TTtot(p) ). The lastcolumn, SseqF F , gives the speedup obtained by comparing against a straight
1 (1)forward sequential FF algorithm (S seqF F (p) = TTtot(p) ).
The results in column K of Table 2 show that, in general, the numberof conflicts that arise in Phase 1 is small and grows as a function of the13δ(p−1)
e
2

0
2
5
11
17

T1
1190
1130
430
260
200

T2
1010
970
280
200
130

T3
0
0
0
0
0

Ttot
2200
2100
710
460
330

Spar
1
1.1
3.1
4.8
6.7

SseqF F
0.6
0.6
1.7
2.6
3.6

0
2
4
0
12

0
2
5
11
17

4400
2250
1300
630
430

3400
1600
1000
800
480

0
0
0
0
0

7800
3850
2300
1430
910

1
2
3.4
5.5
8.6

0.6
1.1
1.9
3.1
4.8

11
12
12
12
12

0
4
12
36
42

0
7
20
46
72

100
55
40
28
20

80
40
20
15
15

0
0
0
0
0

180
95
60
43
35

1
2
3
4.2
5.2

0.6
1.1
1.7
2.3
2.9

13
13
14
13
13

13
13
14
13
13

0
2
14
16
36

0
8
23
53
83

200
130
80
48
40

180
120
50
26
20

0
0
0
0
0

380
250
130
74
60

1
1.5
3
5
6.4

0.5
0.8
1.5
2.7
3.3

1
2
4
8
12

122
142
137
129
121

122
142
137
129
124

0
30
94
94
78

0
155
464
1082
1700

200
110
69
53
55

290
140
72
44
90

0
0
0
1
1

490
250
141
97
145

1
2
3.5
5.6
3.4

0.4
0.8
1.4
2.1
1.4

1
2
4
8
12

377
382
400
407
399

377
382
400
407
399

0
68
98
254
210

0
553
1659
3871
6083

9200
5160
2600
1590
1090

13200
8040
4080
2280
1420

0
3
4
11
8

22400
13203
6684
3881
2518

1
1.7
3.4
5.8
9

0.4
0.7
1.4
2.4
3.7

Problem
mrng2
mrng2
mrng2
mrng2
mrng2

p
1
2
4
8
12

χ1
5
5
5
5
5

χ3
5
5
5
5
5

K
0
0
0
8
18

mrng3
mrng3
mrng3
mrng3
mrng3

1
2
4
8
12

5
5
5
5
5

5
5
5
5
5

598a
598a
598a
598a
598a

1
2
4
8
12

11
12
12
12
12

m14b
m14b
m14b
m14b
m14b

1
2
4
8
12

dense1dense1dense1dense1dense1dense2dense2dense2dense2dense2dTable 2: Experimental results for Algorithm 1.
number of blocks (or processors) p. This agrees well with the result from
Lemma 3.1. We see that for the relatively dense graphs the actual numberof conflicts is much less than the bound given by Lemma 3.1.
The run times obtained show that Algorithm 1 performs as predictedby Lemma 3.2. Particularly, the time required for re-coloring incorrectly
colored vertices is observed to be practically zero (in the order of a few
microseconds) for all our test graphs. This is not surprising as the obtained
value of K is negligible compared to the number of vertices in a given graph.
As results in columns T1 and T2 indicate, the time used to detect conflicts
is approximately the same as the time used to do the initial coloring. This
makes the running time of the algorithm using one processor approximately
double that of the sequential FF. This in turn reduces the speedup obtained
compared to the sequential FF by a factor of 2. The speedup obtained
compared to running the parallel algorithm on one processor gets its best
14

Problem
mrng2
mrng2
mrng2
mrng2
mrng2

p
1
2
4
8
12

χ1
5
5
5
5
5

χ2
5
5
5
5
5

χ4
5
5
5
5
5

K1
0
0
2
16
12

K2
0
0
0
0
0

T1
1050
950
470
300
250

T2
1700
1350
840
500
400

T3
820
650
310
200
170

T4
0
0
0
0
0

Ttot
3570
2650
1620
1000
820

Spar
1
1.4
2.2
3.6
4.4

S2seqF F
0.8
1.0
1.7
2.8
3.4

mrng3
mrng3
mrng3
mrng3
mrng3

1
2
4
8
12

5
5
5
5
5

5
5
5
5
5

5
5
5
5
5

0
0
0
4
24

0
0
0
0
0

3700
1890
1100
540
450

9500
4100
2700
1800
1900

2600
1200
750
450
300

0
0
0
0
0

15800
7190
4550
2790
2650

1
2.2
3.5
5.6
6

0.8
1.8
2.9
4.7
5.0

598a
598a
598a
598a
598a

1
2
4
8
12

11
12
11
12
12

10
10
10
11
11

10
10
10
11
11

0
14
22
40
50

0
0
0
0
0

100
65
35
30
30

200
105
90
99
110

75
37
20
25
15

0
0
0
0
0

375
207
145
154
155

1
1.8
2.6
2.4
2.4

0.8
1.5
2.1
2.0
2.0

m14b
m14b
m14b
m14b
m14b

1
2
4
8
12

13
13
14
13
13

11
12
12
12
11

11
12
12
12
11

0
2
6
12
22

0
0
0
0
0

200
105
70
45
53

520
240
160
120
150

190
80
40
25
20

0
0
0
0
0

910
425
270
190
223

1
2.1
3.4
4.8
4

0.8
1.7
2.7
3.8
3.2

dense1dense1dense1dense1dense1

1
2
4
8
12

122
135
132
126
123

122
122
122
122
121

122
122
122
122
122

0
26
40
104
150

0
0
0
0
2

180
100
80
70
40

250
180
100
80
760

180
140
70
30
30

0
0
0
0
0

610
420
250
180
830

1
1.5
2.5
3.4
0.7

0.7
1.0
1.7
2.4
0.5

dense2dense2dense2dense2dense2

1
2
4
8
12

377
376
394
398
399

376
376
376
376
376

376
376
376
376
376

0
66
112
164
232

0
0
0
0
2

9920
5200
2700
2000
1100

13700
6220
3600
2000
1700

7500
4200
2100
1800
900

0
0
0
0
0

31120
15620
8400
5800
3700

1
2
3.7
5.4
8.4

0.8
1.5
2.8
4.0
6.4

Table 3: Experimental results for Algorithm 2.
values for the two largest graphs mrng3 and dense2.
We see that the number of colors used by Algorithm 1 varies with thenumber of processors used. Comparing column χ 3 of Table 2 and column
χF F of Table 1, we see that the deviation of χ 3 from χF F is at most 1 forgraphs from Problem Set I and II and at most 16% for the two graphs from
Problem Set III.
Algorithm 2. Table 3 lists results of the asynchronous version of Algorithm 2. The number of colors used at the end of Phases 1 and 2 are listed
in columns χ1 and χ2 , respectively. The coloring at the end of Phase 2 is
not guaranteed to be conflict-free. Phases 3 and 4 detect and resolve any
remaining conflicts. Column χ4 lists the number of colors used at the endof Phase 4. The number of conflicts at the end of Phases 1 and 2 are listed
under K1 and K2 , respectively. The time elapsed (in milliseconds) in the15various stages are given in columns T 1 , T2 , T3 , T4 , and Ttot . In order not to
obscure speedup results, the time required to build the color classes prior
to Phase 2 is not included in T2 . In general the time used for this purpose
is in the order of 20% of T1 . Speedup values in column Spar are calculated as in the corresponding column of Table 2. The column S 2seqF F gives
speedups as compared to Culberson’s IG restricted only to two iterations
2 (1)(S2seqF F = T1 (1)+T
Ttot (p) ).
Results in column χ2 confirm that Phase 2 of Algorithm 2 reduces thenumber of colors used by Phase 1. This is especially true for test graphs from
Problem Sets II and III, which contain relatively denser graphs than Problem
Set I. Comparing the results in column χ 2 with the results in columns χF F
and χIDO of Table 1, we see that in general the quality of the coloringobtained using Algorithm 2 is never worse than that of sequential FF and in
most cases comparable with that of the IDO algorithm. IDO is known to be
one of the most effective coloring heuristics [4]. We also observe that, unlike
Algorithm 1, the number of colors used by Algorithm 2 remains reasonably
stable as the number of processors is increased.
From column K2 we see that the number of conflicts that remain after
Phase 2 of Algorithm 2 is zero for almost all test graphs and values of p.
The only occasion where we obtained a value other than zero for K 2 wasusing p = 12 for the graphs dense1 and dense2. These results agree wellwith the claim in Lemma 3.4.
Figure 5 shows the speedup obtained for the problem dense2 using Algorithm 1 and 2 and how the obtained results compare with the ideal speedups.

6

Conclusion

We have presented two new parallel graph coloring heuristics suitable forshared memory programming and analyzed their performance using the
PRAM model.
The heuristics have been implemented using OpenMP and experiments
conducted on an SGI Cray Origin 2000 super computer using very large
graphs validate the theoretical analysis.
The first heuristic is fast, simple and yields reasonably good speedup forgraphs of practical interest run on a realistic number of processors. Generally, the number of colors used by this heuristic never exceeds ∆ + 1. For
relatively dense graphs, the number of colors used by the heuristic increases
slightly as more processors are applied.
The second heuristic is relatively slower, yields reasonable speedup andimproves on the quality of coloring obtained from the first one in that it uses
16

12
s=p
s=p/2
s=2p/3spar(Alg1)
sseqFF(Alg1)
spar(Alg2)
s2seqFF(Alg2)

10speedup (s)

8

6

4

2

0

0

2

4
6
8
number of processors (p)1012Figure 5: Comparison of speedup curves for dense2.
fewer colors. The number of colors used is also more stable as the numberof processors is increased. For the test graphs used in this experiment, thenumber of colors used by this heuristic is in most cases comparable withthat of sequential IDO.
One of the main arguments against using OpenMP has been that it does
not give as good speedup as a more dedicated message passing implementation using MPI. The results in this paper show an example where theopposite is true, the OpenMP algorithms have better speedup than existing
message passing based algorithms. Moreover, implementing the presented
algorithms in a message passing environment would have required a considerable effort and it is not clear if this would have led to efficient algorithms.
Implementing these algorithms using OpenMP is a relatively straight forward task as all the communication is hidden from the programmer.
We point out that in a recent development the algorithms presented inthis paper have been adapted to the CGM model [9].
We believe that the general idea in these coloring heuristics of allowing

17

inconsistency for the sake of concurrency can be applied to develop parallel
algorithms for other graph problems and we are currently investigating thisin problems related to sparse matrix computations.
Acknowledgements We thank the referees for their helpful comments
and George Karypis for making the test matrices in Problem Sets I and II
available.

References[1] J.R. Allwright, R. Bordawekar, P.D. Coddington, K. Dincer, and C. L.
Martin. A comparison of parallel graph coloring algorithms. Technical
Report SCCS-666, Northeast Parallel Architecture Center, Syracuse
University, 1995.
[2] D. Brelaz. New methods to color the vertices of a graph. Communications of the ACM, 22(4), 1979.
[3] G.J. Chaitin, M. Auslander, A.K. Chandra, J. Cocke, M.E. Hopkins,
and P. Markstein. Register allocation via coloring. Computer Languages, 6:47–57, 1981.
[4] T.F. Coleman and J.J. Moré. Estimation of sparse Jacobian matrices
and graph coloring problems. SIAM Journal on Numerical Analysis,
20(1):187–209, 1983.
[5] J.C. Culberson. Iterated greedy graph coloring and the difficulty landscape. Technical Report TR 92-07, Department of Computing Science,
The University of Alberta, Edmonton, Alberta, Canada, June 1992.
[6] A. Gamst. Some lower bounds for a class of frequency assignment problems. IEEE Transactions of Vehicular Technology, 35(1):8–14, 1986.
[7] M.R. Garey and D.S. Johnson. Computers and Intractability. W.H.
Freeman, New York, 1979.
[8] M.R. Garey, D.S. Johnson, and H.C. So. An application of graph coloring to printed circuit testing. IEEE Transactions on Circuits andSystems, 23:591–599, 1976.
[9] A.H. Gebremedhin, I.G. Lassous, J. Gustedt, and J.A. Telle. Graph
coloring on a coarse grained multiprocessor. To be presented at WG
2000, 26th International Workshop on Graph-Theoretic Concepts in
Computer Science, June 15–17, 2000, Konstantz, Germany.
18

[10] G.R. Grimmet and C.J.H. McDiarmid. On coloring random graphs.
Mathematical Proceedings of the Cambridge Philosophical Society,
77:313–324, 1975.
[11] M.T. Jones and P.E. Plassmann. A parallel graph coloring heuristic.
SIAM Journal of Scientific Computing, 14(3):654–669, May 1993.
[12] R.K. Gjertsen Jr., M.T. Jones, and P. Plassman. Parallel heuristics for
improved, balanced graph colorings. Journal of Parallel and Distributed
Computing., 37:171–186, 1996.
[13] G. Karypis. Private Communication.
[14] G. Lewandowski. Practical Implementations and Applications Of Graph
Coloring. PhD thesis, University of Wisconsin-Madison, August 1994.
[15] M. Luby. A simple parallel algorithm for the maximal independent set
problem. SIAM Journal on Computing, 15(4):1036–1053, 1986.
[16] F. Manne. A parallel algorithm for computing the extremal eigenvaluesof very large sparse matrices (extended abstract). Lecture Notes in
Computer Science, Springer, 1541:332–336, 1998.
[17] OpenMP. A proposed industry standard API for shared memory programming. http://www.openmp.org/.
[18] D.J.A. Welsh and M.B. Powell. An upper bound for the chromatic number of a graph and its application to timetabling problems. Computer
Journal, 10:85–86, 1967.

19

Graph Coloring in Optimization Revisited
Assefaw Hadish Gebremedhin ∗Fredrik Manne

Alex Pothen†Abstract
We revisit the role of graph coloring in modeling a variety of matrix partitioning problems that arise in numerical determination oflarge sparse Jacobian and Hessian matrices. The problems consideredin this paper correspond to the various scenarios under which a matrix
computation, or estimation, may be carried out, i.e., the particular
problem depends on whether the matrix to be computed is symmetric or nonsymmetric, whether a one-dimensional or a two-dimensional
partition is to be used, whether a direct or a substitution based evaluation scheme is to be employed, and whether all nonzero entries ofthe matrix or only a subset need to be computed. The resulting complex partitioning problems are studied within a unified graph theoretic
framework where each problem is formulated as a variant of a coloring
problem. Our study integrates existing coloring formulations with new
ones. As far as we know, the estimation of a subset of the nonzero entries of a matrix is investigated for the first time. The insight gained
from the unified graph theoretic treatment is used to develop and analyze several new heuristic algorithms.
Key words: Sparsity, symmetry, Jacobians, Hessians, finite differences, automatic differentiation, matrix partitioning problems, graphcoloring problems, NP-completeness, approximation algorithms∗

Department of Informatics, University of Bergen, N-5020 Bergen, Norway. {assefaw,
fredrikm}@ii.uib.no†
Department of Computer Science, Old Dominion University, Norfolk, VA 23529-0162
USA. pothen@cs.odu.edu and ICASE, NASA Langley Research Center, Hampton, VA
23681-2199 USA. pothen@icase.edu

11Introduction

Algorithms for solving nonlinear systems of equations and numerical optimization problems that rely on derivative information require the repeatedestimation of Jacobian and Hessian matrices. Since this usually constitutes
an expensive part of the entire computation, efficient methods for estimating these matrices via finite difference (FD) or automatic differentiation(AD) techniques are needed. It should be noted that FD techniques find
an approximation whereas AD techniques enable exact (within the limits ofmachine precision) computation. For brevity, we use the term ‘estimation’in referring to both cases. In applying AD and FD techniques, if the sparsity structure of the desired matrix is known a priori, or can be determinedeasily, the nonzero entries can be estimated efficiently. The objective in suchan efficient estimation is to minimize the number of function evaluations or
AD passes required.
This objective calls for a variety of matrix partitioning problems. The
particular problem depends on whether the required matrix is symmetric ornonsymmetric, whether a one-dimensional partition (involving only columns
or rows) or a two-dimensional partition (involving both columns and rows)
is used, whether the entries are evaluated using a direct method or via
substitution, and finally, whether all the nonzero entries of the matrix oronly a subset of them need to be determined.
Several studies have demonstrated the usefulness of graph coloring inmodeling the matrix partitioning problems of our concern [4, 7, 8, 9, 17, 18,23]. However, these studies have been rather disintegrated: a typical study
in this field focuses on one type of matrix, a specific numerical method, and aparticular evaluation scheme. This has at least two consequences. First, theinherent similarity among the various partitioning problems gets obscured.
Second, it makes the identification of a generic formulation difficult thereby
hindering the development of algorithms and software in a flexible manner.
The main purpose of this paper is to study the various matrix partitioning problems within a unified graph theoretic framework. We consider
eight different partitioning problems. For each matrix partitioning problem,
we develop an equivalent graph coloring formulation. To our knowledge,
the problems in partial matrix estimation, where a specified subset of thenonzero entries is to be determined, are studied for the first time. In fullmatrix estimation, the case in which all the nonzero entries are to be determined, we use known coloring formulations for all problems except for theestimation of a nonsymmetric matrix using a one-dimensional partition via
a direct method. For this problem, we propose distance-2 graph coloring as2

Matrix
Jacobian
HessianJacobian
Hessian

1D Partition
distance-2 coloring
distance- 23 coloringNA
acyclic coloring

2D Partition
distance- 23 bicoloringNA
acyclic bicoloringNA

MethodDirectDirectSubstitution
Substitution

Table 1: Graph coloring formulations for estimating all nonzero entries ofderivative matrices. The Jacobian and the Hessian are represented by theirbipartite and adjacency graphs, respectively. NA stands for not applicable.
an alternative formulation to the known distance-1 coloring formulation [7].
The motivation for the distinction between full and partial matrix estimation is the fact that the computation in the latter case can be carried out
even more efficiently.
Table 1 summarizes the coloring formulations used in this paper. The
formulations for one-dimensional estimation of the Hessian are due to Coleman and Moré [8] and Coleman and Cai [4] and those for two-dimensional
estimation of the Jacobian are due to Coleman and Verma [9].
In our work, for all the matrix partitioning problems, we rely on a bipartite graph representation for a nonsymmetric matrix and an adjacency graphrepresentation for a symmetric matrix. We show that these representations
are robust and flexible—they are decoupled from the eventual technique tobe employed and the matrix entries to be computed. This is in contrast withthe column intersection graph representation [7] of a nonsymmetric matrix
which targets at determining all the nonzero entries of the matrix using a
one-dimensional column partition. Moreover, the space required for storing
the column intersection graph of matrix A is proportional to the numberof non-zeros in AT A, whereas the space required for storing the bipartite
graph of A is proportional to the number of non-zeros in A.
Using our graph representations, we identify the distance-2 graph coloring problem as a unifying generic model for the various one-dimensional
matrix partitioning problems, i.e., the coloring problem in a particular caseis some relaxed variant of the distance-2 coloring problem. In the case of
two-dimensional partitioning problems, we build upon the known relationship to graph bicoloring and observe a connection to finding a vertex coverin a graph.
All of the problems listed in Table 1 are known to be NP-hard [4, 8, 9, 21].
We use the insight gained from the interrelationship among these problems
to develop several simple heuristic algorithms for finding sub-optimal solutions.

3

The rest of this paper is organized as follows. Section 2 is a detailed introduction to the partitioning problems that arise in full matrix estimation. In
Section 3 we rigorously develop the equivalent graph problem formulations
and discuss their interrelationship. Section 4 deals with the graph theoretic
formulation of partitioning problems in partial matrix estimation. In Section 5 we compare formulations based on bipartite graph with formulations
based on column intersection graph. In Section 6 we present and analyze
various greedy heuristic algorithms. We conclude the paper in Section 7
with some remarks and point out avenues for further work.

2

Background

2.1

Finite Differences and Partitioning Problems

Given a continuously differentiable function F : R n → Rm , the Jacobian of
F at the point x is the m × n matrix whose (i, j) entry J(x) ij = F 0 (x)ij =∂fi∂xj (x), where f1 (x), f2 (x), . . . , fm (x) are the components of F (x). Let A
denote the Jacobian matrix F 0 (x). An estimate for the jth column of A,
denoted henceforth by aj , can be obtained from the forward difference approximation,
Aej = aj =1∂
F (x) ≈ [F (x + hej ) − F (x)],∂xjh1 ≤ j ≤ n,(1)

where ej is the jth unit vector and h is a positive step length. Other finite
difference approximations of higher order, such as central differences, could
also be used to estimate A. In any case, if F (x) is already evaluated, anapproximation to aj is obtained with one additional function evaluation.
Thus, if each column of A is computed independently, n additional function
evaluations will be required. However, by exploiting the sparsity structure ofA, the required number of function evaluations can be reduced significantly.
The sparsity structure of A is often easily available and the goal here is
to exploit this to estimate the nonzero entries of A using as few function
evaluations as possible under the assumption that evaluating F (x) is more
efficient than evaluating the components f i (x), 1 ≤ i ≤ m, separately.
Given a twice continuously differentiable function f : R n → R, theHessian of f at the point x is the n × n symmetric matrix whose (i, j)∂ 2 f (x). When ∇f is available, ∇2 f can beentry H(x)ij = ∇2 f (x)ij = ∂xi ∂xj
approximated by applying Formula (1) to the function F = ∇f . Again,
we assume that evaluating the gradient ∇f (x) as a single entity is moredesirable than evaluating the components ∂ 1 f (x), . . . , ∂n f (x) separately.
4Let di be the binary vector obtained by adding some unit vectors e j ,j ∈ {1, 2, . . . , n}, together. The problem of estimating a sparse derivative
matrix, Jacobian or Hessian, using FD can then be stated as follows. Given
the sparsity structure of a matrix A find vectors d 1 , d2 , . . . , dp such thatthe products Ad1 , Ad2 , . . . , Adp enable the determination of all the nonzero
entries of A.
Specifying the vectors Ad1 , Ad2 , . . . , Adp gives rise to a system of linear
equations where the unknowns are the nonzero elements of A. If the choiceof the vectors di is such that the resulting system of equations can be ordered
to a diagonal form, then we say that A is directly determined by the vectors
di . If, on the other hand, the vectors d i are chosen such that the system of
equations can be ordered to a triangular form, then the unknowns can be
determined via substitution.
In both a direct and a substitution based determination, minimizing thenumber of function evaluations corresponds to minimizing the number ofvectors p. There is a trade-off in the choice of methods. A direct method ismore restrictive and hence requires more function evaluations compared toa substitution method. On the other hand, a substitution method is subject
to numerical instability, whereas a direct method is not. Moreover, in termsof parallel computation, direct methods offer straightforward parallelization
since the estimates can be read off directly from each row of a matrix-vector
product, whereas substitution methods have less parallelism since there are
more dependencies among the computations required to obtain the matrix
entries. In the next two paragraphs (Sections 2.1.1 and 2.1.2), we consider
minimizing p in a direct and substitution based evaluation, respectively.
2.1.1

Direct estimation

In a direct determination of a matrix A, note that for each nonzero element
aij , there is a vector d in the set {d1 , d2 , . . . , dp } such that aij = (Ad)i , where(Ad)i is the ith component of the vector Ad. Thus, each nonzero matrix
element aij can be read off from some component of the vector Ad.
The problems that arise in the direct, efficient estimation of sparse Jacobian and Hessian matrices can thus be stated as follows.
Problem 2.1 Given the sparsity structure of a general matrix A ∈ R m×n ,find the fewest vectors d1 , d2 , . . . , dp such that Ad1 , Ad2 , . . . , Adp determine
A directly.
Problem 2.2 Given the sparsity structure of a symmetric matrix A ∈ R n×n ,find the fewest vectors d1 , d2 , . . . , dp such that Ad1 , Ad2 , . . . , Adp determine
A directly.
5Curtis et al. [12] were the first to observe that, while using a direct
method, a group of columns can be determined by one evaluation of Ad if
no two columns in this group have a nonzero in the same row position. Suchcolumns are structurally orthogonal, since their pairwise inner products are
zero. Powell and Toint [26] later showed that, in the case of Hessian estimation, the number of function evaluations can be reduced further by considering symmetry. The following two notions define the underlying partitions
used in these methods [8].
Definition 2.3 A partition of the columns of a matrix A is said to beconsistent with a direct determination of A if whenever a ij is a nonzero
element of A then the group containing a j has no other column with anonzero in row i.
Definition 2.4 A partition of the columns of a symmetric matrix A is symmetrically consistent with a direct determination of A if whenever a ij is anonzero element of A then either the group containing a j has no other column with a nonzero in row i, or the group containing a i has no other column
with a nonzero in row j.
Let {C1 , C2 , . . . , Cp } be a consistent partition. With each group C k ,associate a binary vector dk having components δj = 1 if aj belongs to Ck ,and δj = 0 otherwise. Then,
Adk =

X

aj .

aj ∈Ck

If aij 6= 0 and column aj ∈ Ck , then aij = (Adk )i . Thus, all the nonzero
entries of A can be determined with p evaluations of Ad k .
When A is symmetric, a symmetrically consistent partition is sufficient
to determine A directly: if aj is the only column in its group with a nonzero
in row i then aij can be determined as discussed above; alternatively, if a iis the only column in its group with a nonzero in row j then a ji can be
determined.
If a consistent partition (rather than a symmetrically consistent one) isused to compute a symmetric matrix A, the estimate for a ij may actually
be different from that of aji due to truncation error. Thus, using a symmetrically consistent partition to compute half of the nonzero elements ofa matrix and determining the other half by symmetry is preferable bothin terms of reducing computational work and ensuring that the computed
matrix is indeed symmetric.
6

Using Definitions 2.3 and 2.4, Problems 2.1 and 2.2 can be restated asfollows. In the remainder of this paper, we shall use the acronym MPP torefer to a matrix partitioning problem.
Problem 2.5 (MPP1 ) Given the sparsity structure of a matrix A ∈ R m×n ,find a consistent partition of the columns of A with the fewest number ofgroups.
Problem 2.6 (MPP2 ) Given the sparsity structure of a symmetric matrix
A ∈ Rn×n , find a symmetrically consistent partition of the columns of A withthe fewest number of groups.
2.1.2

Estimation via substitution

In estimating a matrix A via a substitution method, the vectors d 1 , . . . , dp
are chosen such that the system of equations defined by the products
Ad1 , . . . , Adp allows the determination of the unknowns via a substitution
process. A partition suitable for a substitution method needs to fulfill morerelaxed requirements compared to a direct method and hence results insmaller number of groups. In the FD context, this fact has been especiallyused when estimating a symmetric matrix since substitution can be effectively combined with the exploitation of symmetry [4]. For a nonsymmetricmatrix, the advantage offered by a substitution method over a direct method
is not so pronounced. An example of a substitution method for a nonsymmetric matrix can be found in [18].
In this paper, we concentrate on using a substitution method for a symmetric matrix. To illustrate the fact that a partition used in a substitution
method requires fewer groups than that used in a direct method, consider
the 4 × 4 symmetric matrix H shown in Figure 1.
X X
X X X
X X X
X X

Figure 1: The structure of a 4 × 4 symmetric matrix
For matrix H, any partition consistent with a direct determination requires at least three groups. One such partition is {(h 1 , h4 ), (h2 ), (h3 )},where hj is the jth column of H. However, if we do not insist on determining the elements directly, two groups would suffice. For example{(h1 , h3 ), (h2 , h4 )}. The two matrix-vector products corresponding to the7two groups yield a system of eight equations involving the nonzero entries of
H. Note that, due to symmetry, nonzero element h ij can be identified withhji , and hence, there are effectively seven unknowns in the system. Thissystem can be ordered to a triangular form and be solved via substitution.
In general, a partition of the columns of a symmetric matrix induces asubstitution method if there is an ordering of the matrix unknowns suchthat all unknowns can be solved for, in that order, using symmetry andpreviously solved elements.
We now formally define such a partition and then state the corresponding
partitioning problem.
Definition 2.7 A partition of the columns of a symmetric matrix A is said
to be substitutable if there exists an ordering on the elements of A such thatfor every nonzero aij , either aj is in a group where all the nonzeros in row
i, from other columns in the same group, are ordered before a ij or ai is ina group where all the nonzeros in row j, from other columns in the samegroup, are ordered before aij .
Problem 2.8 (MPP3 ) Given the sparsity structure of a symmetric matrix
A ∈ Rn×n , find a substitutable partition of the columns of A with the fewest
number of groups.

2.2

Automatic Differentiation and Partitioning Problems

Automatic differentiation (AD) is a chain rule based technique for evaluating
the derivatives of functions defined by computer programs. The two basic
modes of operation of AD, known as forward and reverse mode, correspondto a bottom-up and a top-down strategy of accumulating partial derivativesof elementary functions that define the computational scheme of the function
to be differentiated. A treatment of the technical details of AD is beyond
the scope of this paper, but the interested reader is referred to, for instance,
the books [14] and [10].
What is of interest for us is that, as in the FD setting, the efficientcomputation of matrices using AD gives rise to partitioning problems inwhich structural orthogonality continues to be the partition-criterion. In
particular, one can use the forward mode to compute a group of columnsof a matrix A from the product Ad1 , where the vector d1 has nonzeros in
positions corresponding to columns of A that are structurally orthogonal.
Furthermore, in the reverse mode, a group of structurally orthogonal rowsof A can be computed from the vector-matrix product d 2 T A, where d2 is8

an appropriately defined vector. This means that one can potentially takeadvantage of the sparsity available in columns and in rows.
One way of exploiting sparsity in columns and rows is to separately partition the columns and rows of the matrix and use the partition which gives
the minimum number of groups. For a symmetric matrix, a row partitionis equivalent to a column partition, but for a nonsymmetric matrix, the twopartitions may differ considerably. For example, consider an n × n matrix
where all the entries on the diagonal and the first row are nonzero, and therest of the matrix entries are all zero. For such a matrix structure, a column
partition requires n groups whereas a row partition requires just two groups.
However, an approach based on a separate row and column partition is
not always satisfactory. For example, consider an n × n matrix where allof the elements in the first row, first column, and the diagonal are nonzero
and the rest of the entries are all zero. For such a structure, a row partition
requires n groups and so does a column partition. However, using a combined
row and column partition, three groups are enough to determine all thenonzero entries of the matrix. First, separately evaluate the entries in thefirst column and the first row (two groups). Then, since the remaining(n − 1) × (n − 1) matrix is diagonal, determine all entries by one evaluation.
Thus, three groups (two column and one row) suffice to determine all thenonzero entries.
In Sections 2.2.1 and 2.2.2 we consider such a computation of a nonsymmetric matrix using the combined modes of AD via direct and substitution
methods. We call a partition that involves both rows and columns a twodimensional partition as opposed to a one-dimensional partition in which
either only rows or only columns are involved. Note that a two-dimensional
partition does not make sense for computing a symmetric matrix. In particular, a symmetry-exploiting one-dimensional partition is sufficient.
2.2.1

Direct computation

Consider the vectors d1 , d2 , . . . , dp in Problem 2.1 as the p columns of then × p matrix D. An alternative way of posing the problem would then be:
given the structure of a matrix A ∈ R m×n , find a matrix D ∈ Rn×p withthe least value of p such that the product AD determines A directly. By thesame token, the problem that arises in the two-dimensional efficient directcomputation of a Jacobian can be posed as follows.
Problem 2.9 Given the sparsity structure of the matrix A ∈ R m×n , findmatrices D1 ∈ Rn×p1 and D2 ∈ Rm×p2 such that AD1 and D2T A together
determine A directly and the value p = p 1 + p2 is minimized.
9Hossain and Steihaug [17] studied Problem 2.9 and reformulated it as apartitioning problem by using the notion of consistent row-column partitionin which the entire set of rows and columns is partitioned into two respective set of groups. Coleman and Verma [9] also studied the same problem
and identified a similar two-dimensional partition problem. Their notion ofpartition differs from that of Hossain and Steihaug in that it partitions only
a subset of the columns and the rows of the matrix that suffice for the direct
determination of the entries. The following concepts were used to formalize
the requirements.
Definition 2.10 A bipartition of a matrix A is a pair (Π C , ΠR ) where ΠC is
a column partition of a subset of the columns of A and Π R is a row partitionof a subset of the rows of A.
Definition 2.11 A bipartition (ΠC , ΠR ) of a matrix A is consistent with adirect determination if for every nonzero entry a ij of A, either column j is ina group of ΠC which has no other column having a nonzero in row i, or row
i is in a group of ΠR which has no other row having a nonzero in column j.

The number of column and row groups in a consistent bipartition corresponds to the number of forward and reverse AD passes, respectively,
required to compute the nonzero entries directly. To see this, observe thata nonzero aij can be determined either from a column group where column
j is the only column with a nonzero in row i, or from a row group whererow i is the only row with a nonzero in column j. Hence, assuming that thecomputational costs involved in the forward and reverse modes of AD areof the same order, in an efficient method, the value |Π C | + |ΠR | is required
to be as small as possible1 .
Thus Problem 2.9 can be reformulated as follows.
Problem 2.12 (MPP4 ) Given the sparsity structure of a matrix A ∈Rm×n , find a bipartition (ΠC , ΠR ) of A consistent with a direct determination such that |ΠC | + |ΠR | is minimized.
2.2.2

Computation via Substitution

In using a substitution method in the AD context, the requirement on the
bipartition can be relaxed so as to obtain fewer number of groups compared1

In this paper, we are concerned only with computational cost; however, in general,
the forward mode requires less memory space than the reverse mode.

10

Matrix
Jacobian
HessianJacobian
Hessian

1D Partition
MPP1
MPP2
NA
MPP3

2D Partition
MPP4
NA
MPP5
NA

MethodDirectDirectSubstitution
Substitution

Table 2: Partitioning problems in estimating/computing all nonzero entriesof derivative matrices using FD and AD. The entry NA denotes that thecase is not applicable.
with that in a direct method. We state the following definition used by
Coleman and Verma [9] to subsequently give the fifth matrix partitioning
problem of our concern.
Definition 2.13 A bipartition (ΠC , ΠR ) of a matrix A is consistent with adetermination by substitution if there exists an ordering on the elements ofA such that for every nonzero entry a ij , either column j is in a group where
all nonzeros in row i, from other columns in the group, are ordered before
aij or row i is in a group where all the nonzeros in column j, from otherrows in the group, are ordered before a ij .
Problem 2.14 (MPP5 ) Given the sparsity structure of a matrix A ∈Rm×n , find a bipartition (ΠC , ΠR ) of A consistent with a determination
by substitution such that |ΠC | + |ΠR | is minimized.
Problems MPP1 through MPP5 are summarized in Table 2. Note thatthe problems in Table 2 are formulated independent of the numerical technique used. For example, MPP1 could arise in using FD or only the forward
mode of AD.

2.3

Other Methods

The matrix estimation methods considered in this paper rely on a onedimensional or a two-dimensional partition (symmetrically) consistent witha direct or a substitution-based determination. However, approaches where
this is not necessarily the case have also been suggested. For instance,
direct methods that allow structurally non-orthogonal columns to reside in
the same group and/or allow columns to reside in several groups have beensuggested [8, 25, 26]. McCormick [23] gives a classification of direct methods
for estimating a symmetric matrix, including those that do not necessarily
rely on consistent partitions.

11In another direction, Hossain [16] suggests a technique for a direct estimation of a Jacobian in which the rows are first grouped into blocks that define ‘segmented’ columns and then the segments are partitioned into groups
each of which consists of structurally independent segments. It is shown
that, for some matrix structures, such an approach may reduce the number of function evaluations compared with an approach that does not use
segmentation.

3

Graph Formulations

Problems MPP1 through MPP5 are combinatorial and several previous studies have demonstrated the usefulness of graph theoretic approaches in analyzing and solving them [4, 7, 8, 9, 17, 23]. In this section, we integratethese approaches within a unified framework. In addition, we propose anew, more flexible, graph formulation for problem MPP1. The graph formulation of each matrix problem depends on the choice of the graph used
to represent the matrix structure. In Section 3.2 we describe our graph representations which will be used to develop the equivalent graph problems
in Sections 3.3 through 3.5. In Section 3.6 we show how the various graphproblems obtained relate with one another. We start this section by definingsome basic graph theoretic concepts. Other graph concepts will be definedlater as required.

3.1

Basic Definitions

A graph G is an ordered pair (V, E) where V is a finite and nonempty
set of vertices and E is a set of unordered pairs of distinct vertices callededges. If (u, v) ∈ E, vertices u and v are said to be adjacent; otherwise
they are called non-adjacent. A path of length l in a graph is a sequence
v1 , v2 , . . . , vl+1 of distinct vertices such that vi is adjacent to vi+1 , for 1 ≤i ≤ l. Two distinct vertices are said to be distance-k neighbors if the shortest
path connecting them has length at most k; otherwise they are called nondistance-k neighbors. The number of distance-k neighbors of a vertex u isreferred to as the degree-k of u.
In a graph G = (V, E), a set of vertices C ⊆ V is said to cover a set ofedges F ⊆ E if for every edge e ∈ F , at least one of the endpoints of e isin C. If the set C covers the entire E, it is called a vertex cover. The setof vertices I ⊆ V is called an independent set if no pair of vertices in I are
adjacent to each other.
A graph G = (V, E) is bipartite if its vertex set V can be partitioned into

12

two disjoint sets V1 and V2 such that every edge in E connects a vertex from
V1 to a vertex from V2 . We denote a bipartite graph by Gb = (V1 , V2 , E).

3.2

Representing Matrix Structures Using Graphs

Bipartite graph Let A be an m × n matrix with rows r 1 , r2 , . . . , rm andcolumns a1 , a2 , . . . , an . We define the bipartite graph Gb (A) of A as Gb (A) =(V1 , V2 , E) where V1 = {r1 , r2 , . . . , rm }, V2 = {a1 , a2 , . . . , an }, and (ri , aj ) ∈E whenever aij is a nonzero element of A, for 1 ≤ i ≤ m, 1 ≤ j ≤ n.
Note that Gb (A) is a space-efficient representation for a nonsymmetricmatrix A. To see this, notice that the number of vertices |V 1 |+|V2 | = m+n,
and the number of edges |E| = nnz(A), where nnz(A) is the number ofnonzeros in A. Also, note that the graph can be constructed by reading off
the entries of the matrix without any further computation.
Adjacency graph Let A ∈ R n×n be a symmetric matrix with nonzero
diagonal elements and let its columns be a 1 , a2 , . . . , an . The adjacency graphof A is G(A) = (V, E) where V = {a1 , a2 , . . . , an }, and (ai , aj ) ∈ E wheneveraij is a nonzero element of A, for 1 ≤ i ≤ n, 1 ≤ j ≤ n, i 6= j. Notethat G(A) is a space-efficient, symmetry-exploiting, graph representationof the symmetric matrix A with no explicit representation for the edges
corresponding to the nonzero diagonal elements. In particular, the numberof vertices is n and the number of edges is 21 (nnz(A)−n). This is in contrast
to 2n and nnz(A), respectively, had a bipartite graph representation beenused.
For the symmetric matrix of our interest, the Hessian matrix, the assumption that the diagonals are nonzero is reasonable in many contexts. In
particular, the Hessian is usually positive definite [4].
To simplify notation, we shall use a i both when referring to the ith
column of the matrix A and the corresponding vertex in G(A) or G b (A).
Column intersection graph Our graph formulations for problems MPP1
through MPP5 are based on the bipartite and adjacency graph representations discussed above. However, in the literature, other graph representations have also been used. In particular, Coleman and Moré [7] used thecolumn intersection graph Gc (A) to represent a nonsymmetric matrix A. In
Gc (A) = (V, E), the columns of A constitute the vertex set V , and an edge
(ai , aj ) is in E whenever columns ai and aj have nonzero entries at the samerow position, i.e. whenever ai and aj are not structurally orthogonal.

13

3.3

Distance-k Graph Coloring

A distance-k p-coloring, or (k, p)-coloring for short, of a graph G = (V, E) isa mapping φ : V → {1, 2, . . . , p} such that φ(u) 6= φ(v) whenever u and v are
distance-k neighbors. The minimum possible value of p in a (k, p)-coloringof a graph G is called its k-chromatic number, and is denoted by χ k (G).
A (k, p)-coloring of G = (V, E) is called partial if it involves only a subsetof the vertices; in particular, a partial (k, p)-coloring of G = (V, E) on W ,W ⊂ V , is a mapping φ : W → {1, 2, . . . , p} such that φ(u) 6= φ(v) wheneveru and v are distance-k neighbors.
The distance-k graph coloring problem (DkGCP) asks for an optimal(k, p)-coloring of a graph: given a graph G and an integer k, find a (k, p)coloring of G such that p is minimized.
Notice that a (k, p)-coloring of G = (V, E) partitions the set V into pgroups (called color classes) U1 , U2 , . . . , Up , where Ui = {u ∈ V : φ(u) = i}.
Each color class is a distance-k independent set, i.e., no pair of distinct
vertices consists of distance-k neighbors. Noting this, a natural question
arises—can the matrix partitioning problems MPP1 through MPP5 be related to the DkGCP for some value of k?
Recall that structural orthogonality is the criterion used in a consistentpartition of the columns (or rows) of a matrix. The following two simple
observations provide the graph theoretic equivalents of structural orthogonality in matrices.
Lemma 3.1 Let A ∈ Rm×n be a matrix and Gb (A) = (V1 , V2 , E) be its
bipartite graph. Two columns (or rows) of A are structurally orthogonal if
and only if the corresponding vertices in G b (A) are non-distance-2 neighbors.

Proof: We prove the statement for columns; a similar argument can beused to prove the case for rows. Assume that vertices a i and aj in V2 are
non-distance-2 neighbors. Thus, by definition, there is no path a i , rk , aj in
Gb for any rk ∈ V1 , 1 ≤ k ≤ m. In terms of matrix A, this means that
there is no k ∈ [1, m] such that both aki and akj are nonzero. Hence, bydefinition, ai and aj are structurally orthogonal.
To prove the ‘only if’ part of the statement, assume that columns a i
and aj are structurally orthogonal. Then, by definition, there is no k ∈[1, m] such that aki 6= 0 and akj 6= 0. This implies that there is no pathai , rk , aj in Gb (A), for any 1 ≤ k ≤ m. Hence, by definition, a i and aj are
non-distance-2 neighbors.
2

14

Similarly, one can prove the following statement for the case of a symmetric matrix and its adjacency graph representation.
Lemma 3.2 Let A ∈ Rn×n be a symmetric matrix with nonzero diagonal
elements and let G(A) = (V, E) be its adjacency graph. Two columns in A
are structurally orthogonal if and only if the corresponding vertices in G(A)are non-distance-2 neighbors.
Lemmas 3.1 and 3.2 provide a partial answer to our question regarding
the relationship between the matrix partitioning problems and the DkGCP.
As discussed in the forthcoming sections, distance-2 coloring is a genericmodel in efficient derivative matrix estimation using methods that rely on
one-dimensional column or row partition.

3.4

Coloring Problems in Direct Methods

In this subsection we consider problems MPP1, MPP2, and MPP4.
By Lemma 3.1, finding a consistent partition of the columns of a matrix
A is equivalent to finding a partial distance-2 coloring of G b (A) = (V1 , V2 , E)
on V2 . The following result formalizes the equivalence.
Theorem 3.3 Let A be a nonsymmetric matrix and G b (A) = (V1 , V2 , E)
be its bipartite graph representation. A mapping φ is a partial distance-2
coloring of Gb (A) on V2 if and only if φ induces a consistent partition of thecolumns of A.
In view of Theorem 3.3, Problem MPP1 is equivalent to the following
graph coloring problem (GCP).
Problem 3.4 (GCP1 ) Given the bipartite graph G b (A) = (V1 , V2 , E) representing the sparsity structure of a matrix A ∈ R m×n , find a partial (2,p)coloring of Gb (A) on V2 with the least value of p.
For matrices with a few dense rows, a row partition may yield fewergroups than a column partition. Consequently, the matrix problem oneneeds to solve is MPP1 applied on AT . In such cases, our graph formulation becomes handy—the equivalent problem is to find a partial distance-2
coloring on the vertex set V1 .
By Lemma 3.2, finding a consistent partition of the columns of a symmetric matrix A is equivalent to finding a distance-2 coloring of the adjacency
graph G(A). This equivalence was in fact first observed by McCormick[23]. However, as has been stated earlier, the symmetry present in A can be15

exploited to further reduce the number of groups (colors) required. Thus,
we now consider the graph coloring formulation of the partitioning problem
where A is symmetric (MPP2).
Consider a symmetric matrix A with nonzero diagonal elements and letaij , i 6= j, be any nonzero element in A. Further, let a ki , k 6= i, j andajl , l 6= i, j, k be any other two nonzero elements. By Definition 2.4, in asymmetrically consistent partition of A,• columns ai and aj should belong to two different groups (this is because
both aii and ajj are nonzero), and• columns aj and ak should belong to two different groups, or columns
ai and al should belong to two different groups.
Coleman and Moré [8] gave an equivalent characterization of the aforementioned conditions in terms of a coloring of the associated adjacency
graph. Specifically, they introduced the notion of distance- 23 coloring 2 , defined below.
Definition 3.5 A mapping φ : V → {1, 2, . . . , p} is a ( 23 , p)-coloring of thegraph G = (V, E) if φ is (1, p)-coloring of G and every path of length threeuses at least three colors.
The name ‘distance- 23 coloring’ is chosen to reflect that it is in a sense‘in-between’ distance-1 and distance-2 colorings. In particular, a distance3
2 coloring is a relaxed distance-2 and a restricted distance-1 coloring. As
an illustration, observe that a distance-1 coloring requires two colors forevery path of length one, a distance-2 coloring requires three colors for every
path of length two, and a distance- 23 coloring is a distance-1 coloring further
restricted to require three colors for every path of length three (see Figure 2).
Note that a 4-cycle requires two, three, and four colors in a distance-1, a
distance- 32 and a distance-2 coloring, respectively.
The following theorem formalizes the connection between symmetricallyconsistent partition and distance- 23 coloring. The result follows directly fromthe discussion that led to the definition of ( 32 , p)-coloring.
Theorem 3.6 [Coleman and Moré [8]]Let A be a symmetric matrix with nonzero diagonal elements and G(A) =(V, E) be its adjacency graph representation. A mapping φ is a ( 23 , p)coloring of G(A) if and only if φ induces a symmetrically consistent partitionof the columns of A.
2

Coleman and Moré used the term path-coloring.

16

1

2

1

2

3

1

2

1

2

2

1

4

3

1

2

1

1

2

2

3

3

Figure 2: Distance-1, 2, and 32 coloring of paths and a 4-cycle.
By Theorem 3.6, the following problem is equivalent to MPP2.
Problem 3.7 (GCP2 ) Given the adjacency graph G(A) = (V, E) representing the sparsity structure of a symmetric matrix A ∈ R n×n with nonzero
diagonal elements, find a ( 23 , p)-coloring of G(A) with the least value of p.
Our next problem, MPP4, aims at finding a bipartition with the fewest
number of groups consistent with a direct determination. When a nonsymmetric matrix A is represented by its bipartite graph G b (A) = (V1 , V2 , E),
we have shown that a column partition consistent with a direct determination can be obtained by finding a partial distance-2 coloring of G b on V2 .
We now consider how this coloring has to be modified to capture a bipartition consistent with a direct determination. Notice that the coloring we arelooking for should meet the following conditions.
• The sets V1 and V2 should use disjoint set of colors.
• Some vertices may not be involved in the determination of any nonzero
entry of the underlying matrix. Such vertices are assigned a ‘neutral’color (say color zero).
• Since every nonzero matrix entry has to be determined, for every edge
in E, at least one of the endpoints has to be assigned a nonzero color.
• A nonzero matrix entry may be determined either from a positively
colored column vertex or a positively colored row vertex. This suggests that the coloring condition sought here is some relaxation of the
distance-2 coloring requirement imposed in the case of one-dimensional
partition.
The following definition, introduced by Coleman and Verma [9], albeit
using a different terminology, formalizes the conditions listed above. The
subsequent theorem establishes the equivalence between the matrix andgraph problems.
17Definition 3.8 Let Gb = (V1 , V2 , E) be a bipartite graph. A mapping
φ : [V1 , V2 ] → {0, 1, . . . , p} is a distance- 32 bicoloring of Gb if the following
conditions hold.
1. If u ∈ V1 and v ∈ V2 , then φ(u) 6= φ(v) or φ(u) = φ(v) = 0.
2. If (u, v) ∈ E, then φ(u) 6= 0 or φ(v) 6= 0.
3. If vertices u and v are adjacent to vertex w with φ(w) = 0, thenφ(u) 6= φ(v).
4. Every path of three edges uses at least three colors.
Theorem 3.9 [Coleman and Verma[9]]Let A be an m×n matrix and Gb (A) = (V1 , V2 , E) be its bipartite graph. The
mapping φ : [V1 , V2 ] → {0, 1, . . . , p} is a distance- 23 p-bicoloring if and only
if φ induces a bipartition (ΠC , ΠR ) of A, with |ΠC | + |ΠR | = p, consistent
with a direct determination.
Thus MPP4 is equivalent to the following graph problem.
Problem 3.10 (GCP4 ) Given the bipartite graph G b (A) = (V1 , V2 , E)
representing the sparsity structure of an m × n matrix A, find a distance- 23
p-bicoloring of Gb (A) with the least value of p.

3.5

Coloring Problems in Substitution Methods

In this subsection, we consider problems MPP3 and MPP5. To formulate
MPP3 as a graph problem, we introduce the notion of acyclic coloring.
Coleman and Cai [4] established the connection between acyclic 3 coloring
and the estimation of a symmetric matrix using a substitution method.
Acyclic coloring had been studied earlier by Grünbaum [15] in a different
context.
Definition 3.11 A mapping φ : V → {1, 2, . . . , p} is an acyclic p-coloringof a graph G = (V, E) if φ is (1, p)-coloring and every cycle in G uses atleast three colors.
Theorem 3.12 [Coleman and Cai [4]]Let A be a symmetric matrix with nonzero diagonal elements and G(A) =(V, E) be its adjacency graph representation. A mapping φ is an acyclic
p-coloring of G(A) if and only if φ induces a substitutable partition of thecolumns of A.
3

Coleman and Cai use the term cyclic coloring to refer to what is known as acyclic
coloring in the graph theoretic literature.

18

Consider the 5 × 5 symmetric matrix and its adjacency graph depicted
in Figure 3. An optimal acyclic coloring that uses three colors is shown. It
can be verified that this induces a substitutable partition of the columns.
By contrast, a symmetric partition consistent with a direct determination(a distance- 32 coloring) would have required four groups (colors).
a

b

c

d

X X

e

2

X

b

X X

X

X

X

X

X

X

X

X

X

X

1
c

1
a

e
2

d
3

Figure 3: A 5 × 5 symmetric matrix and an acyclic coloring of its adj. graphIn view of Theorem 3.12, MPP3 is equivalent to the following graphproblem.
Problem 3.13 (GCP3 ) Given the adjacency graph G(A) = (V, E) representing the sparsity structure of a symmetric matrix A ∈ R n×n with nonzero
diagonal elements, find an acyclic p-coloring of G(A) with the least value ofp.
The relationship between bicoloring and bipartition, established by Theorem 3.9, coupled with that between acyclic coloring and substitutable partition, established by Theorem 3.12, suggests that ‘acyclic bicoloring’ might
be the right graph model for MPP5. Coleman and Verma [9] showed that
this was indeed the case.
Definition 3.14 Let Gb = (V1 , V2 , E) be a bipartite graph. A mapping
φ : [V1 , V2 ] → {0, 1, . . . , p} is an acyclic bicoloring of G b if the following
conditions hold.
1. If u ∈ V1 and v ∈ V2 , then φ(u) 6= φ(v) or φ(u) = φ(v) = 0.
2. If (u, v) ∈ E, then φ(u) 6= 0 or φ(v) 6= 0.
3. If vertices u and v are adjacent to vertex w with φ(w) = 0, thenφ(u) 6= φ(v).
4. Every cycle uses at least three colors.

19

Matrix
Jacobian
HessianJacobian
Hessian

1D Partition
distance-2 coloring
distance- 23 coloringNA
acyclic coloring

2D Partition
distance- 23 bicoloringNA
acyclic bicoloringNA

MethodDirectDirectSubstitution
Substitution

Table 3: Graph coloring formulations for estimating all nonzero entries ofderivative matrices. The Jacobian and the Hessian are represented by theirbipartite and adjacency graphs, respectively. NA stands for not applicable.
Theorem 3.15 [Coleman and Verma [9]]Let A be an m×n matrix and Gb (A) = (V1 , V2 , E) be its bipartite graph. The
mapping φ : [V1 , V2 ] → {0, 1, . . . , p} is an acyclic p-bicoloring if and only if
φ induces a bipartition (ΠC , ΠR ) of A, with |ΠC | + |ΠR | = p, consistent withdetermination by substitution.
By Theorem 3.15, the following coloring problem is equivalent to MPP5.
Problem 3.16 (GCP5 ) Given the bipartite graph G b (A) = (V1 , V2 , E)
representing the sparsity structure of an m × n matrix A, find an acyclic
p-bicoloring of Gb (A) with the least value of p.
Table 3 (the same table was given in Section 1) summarizes the graphcoloring problems that arise in efficient derivative matrix estimation.

3.6

Distance-k Chromatic Numbers

In this subsection, we expose the inter-relationship among the various colorings introduced thus far and show that distance-2 coloring is the mostgeneral among them.
The power of a graph gives an alternative view to the DkGCP. The kth
power of a graph G = (V, E) is the graph G k = (V, F ) where (u, v) ∈ F if
and only if u and v are distance-k neighbors in G. The following equivalence
follows immediately.
Lemma 3.17 Let Gk be the kth power of graph G. A mapping φ is a (k, p)coloring of G if and only if it is a (1, p)-coloring of G k .
A particular implication of Lemma 3.17 is that distance-2 coloring of agraph is equivalent to distance-1 coloring of the square of the graph. Thisestablishes the equivalence between GCP1, our bipartite graph based formulation of MPP1, and the distance-1 coloring formulation, which uses the20column intersection graph, suggested by Coleman and Moré [7]. Specifically,as has been shown in [7], the column intersection graph G c (A) of a matrix
A is isomorphic to the adjacency graph of A T A. We note that Gc (A) is
in fact the subgraph of Gb (A)2 induced by the vertices in V2 . For a graph
G = (V, E), let the graph induced by U ⊆ V be denoted by G[U ].
Lemma 3.18 Let Gb (A) = (V1 , V2 , E) and Gc (A) = (V2 , E 0 ) be the bipartiteand column intersection graphs of matrix A. Then, G c = Gb 2 [V2 ].
Observe that the distance-1 neighbors of a vertex in a graph G form aclique in the square of the graph. A clique is a set of vertices in which thevertices are mutually adjacent to each other. This observation immediately
provides a lower bound on χ2 (G), the 2-chromatic number of G. Let ∆denote the maximum degree-1 in G.
Lemma 3.19 For any graph G, χ2 (G) ≥ ∆ + 1.
Proof: Observe that the cardinality of a maximum clique in the squaregraph G2 is ∆ + 1.
2
Further, the conditions required by distance-1 coloring, acyclic coloring,
distance- 32 coloring, acyclic bicoloring, distance- 23 bicoloring, and distance-2
coloring imply the following relationships among their respective chromatic
numbers. The chromatic numbers for distance-k coloring (and bicoloring)of a general graph G (and bipartite graph G b ) are denoted by χk (G) (andχkb (Gb )). Similarly the chromatic numbers for acyclic coloring (and bicoloring) of a graph G (and a bipartite graph G b ) are denoted by χa (G) (andχab (Gb )).
Theorem 3.20 For a general graph G = (V, E),
χ1 (G) ≤ χa (G) ≤ χ 3 (G) ≤ χ2 (G) = χ1 (G2 ).
2

Proof: Observe that a distance-2 coloring is a distance- 23 coloring; a distance3
2 coloring is an acyclic coloring; and an acyclic coloring is a distance-1
coloring.
2
Theorem 3.21 For a bipartite graph G b = (V1 , V2 , E),
χab (Gb ) ≤ χ 3 b (Gb ) ≤ min{χ1 (Gb 2 [V1 ]), χ1 (Gb 2 [V2 ])},
2

where (Gb 2 [W ]) is the sub-graph of Gb 2 induced by W .
21

Proof: The first inequality is obvious. For the second inequality, observethat a partial distance-2 coloring on V 2 is a valid distance- 23 bicoloring where
all the vertices in V1 are specified to be colored with 0. A similar argument,with the roles of V1 and V2 interchanged, can be used to complete the proof.
2
In the context of matrix estimation using numerical methods, the implication of Theorem 3.21 is that, an optimal two-dimensional partition,
irrespective of the structure of the matrix, yields fewer (or at most as many)groups compared to an optimal one-dimensional partition, and hence potentially results in a more efficient computation.
The results in this section show that distance-2 coloring is an archetypal
model in the estimation of Jacobian and Hessian matrices using techniquesthat rely on one-dimensional partition via direct and substitution methods.
Distance-2 coloring has other applications. Examples include channel
assignment [20] and facility location problems (see Chapter 5 in [27]). Froma more theoretical perspective, distance-2 coloring for planar graphs has
been studied in [2] and a similar study for chordal graphs is available in [1].

4Partial Matrix Estimation

In many PDE constrained optimization contexts, the Jacobian or the Hessian is formed only for preconditioning purposes. For preconditioning, it isoften common to compute a subset of the matrix elements. Computing agood preconditioner is critical for fast convergence to the solution. The recent survey article [19] discusses various applications where methods known
as “Jacobian-free Newton-Krylov” are used. A basic ingredient in the useof these methods is an approximate computation of some elements of theJacobian. Also, there are examples in which only certain elements of the
Hessian need to be updated in an iterative procedure, while others do notbecause they are unlikely to change in value [3].
In this section we develop the graph coloring formulations of the partitioning problems that arise in the estimation of a specified subset of thenonzero entries of a matrix via direct methods. We call this partial matrix
estimation as opposed to full matrix estimation, where all nonzero entries
are required to be determined.
The coloring formulations in this section are new and more sophisticatedthan the coloring formulations in full matrix estimation. The motivation
for developing the new graph formulations is that efficient partial matrix
estimation can be used to further reduce the number of colors needed to22estimate the required elements. For example, if only the diagonal elementsof a Hessian are needed, then we need only the distance-1 coloring of theadjacency graph, rather than the distance- 23 coloring required for full matrix
estimation.
The colorings defined in this section allow a vertex to have the color
zero. A vertex with color zero signifies the fact that it would not be used toestimate any element of the matrix represented by the graph, i.e., columns orrows that correspond to the color zero are not used to estimate any elementsin those columns or rows.
The rest of this section is organized in three parts. Each part deals witha scenario defined by the kind of matrix under consideration (symmetricor nonsymmetric) and the type of partition employed (one-dimensional or
two-dimensional). In each case, the required entries are assumed to be determined using a direct method. The problems that correspond to estimation
via substitution are not considered in this paper.

4.1

Nonsymmetric matrix, One-dimensional partition

Let A ∈ Rm×n be a nonsymmetric matrix, and S denote the set of nonzero
elements of A required to be estimated. A partition {C 1 , . . . , Cp } of a subsetof the columns of A is consistent with a direct determination of S if for every
aij ∈ S, column aj is included in some group that contains no other column
with a nonzero in row i.
Let Gb (A) = (V1 , V2 , E) be the bipartite graph of A, and F ⊆ E correspond to the elements of S. A mapping φ : V 2 → {0, 1, . . . , p} is a distance2 coloring of Gb restricted to F if the following conditions hold for every(v, w) ∈ F , where v ∈ V1 , w ∈ V2 .
1. φ(w) 6= 0, and
2. for every path (u, v, w), φ(u) 6= φ(w).
Theorem 4.1 The mapping φ is a distance-2 p-coloring of G b (A) restricted
to F if and only if φ induces a column partition {C 1 , . . . , Cp } consistent witha direct determination of S.
Proof: Assume that φ is a distance-2 p-coloring of G b (A) restricted to F .
We show that the groups {C1 , . . . , Cp } where Cα = {aj : φ(aj ) = α}, 1 ≤α ≤ p, constitute a partition of a subset of the columns of A consistent witha direct determination of S. First, observe that by coloring condition 1, forevery aij ∈ S ((ri , aj ) ∈ F ), φ(aj ) 6= 0 and thus column aj belongs to groupCφ(aj ) and hence is involved in the partition. Assume now that the partition23

induced by the coloring is not consistent with a direct determination of S.
This occurs only if there exist nonzero elements a ij and aik , k 6= j, such thataij ∈ S and both aj and ak belong to group Cα0 for some α0 , 1 ≤ α0 ≤ p.
But this contradicts coloring condition 2, and hence cannot occur.
Conversely, assume that the partition C = {C 1 , . . . , Cp } is consistent
with a direct determination of S. Construct a coloring φ of G b (A) as follows.
φ(aj ) = α if aj ∈ Cα , and φ(aj ) = 0 if aj 6∈ C. We claim that φ is a distance2 p-coloring of Gb (A) restricted to F . Each vertex in V 2 incident to an edge
in F corresponds to a column with an entry in S and thus gets a nonzero
color. Thus φ satisfies coloring condition 1. Consider any path (a j , ri , ak )where (ri , aj ) ∈ F . Note that such a path in Gb (A) exists whenever entries
aij and aik are nonzero. The partition condition implies that column a kcannot be in the same group as aj . Thus, by construction, φ(aj ) 6= φ(ak ),satisfying coloring condition 2.
2

4.2

Symmetric matrix, One-dimensional partition

Let A ∈ Rn×n be a symmetric matrix with nonzero diagonal elements, and S
denote the set of nonzero elements of A required to be estimated. A partition{C1 , . . . , Cp } of a subset of the columns of A is symmetrically consistent witha direct determination of S if for every a ij ∈ S at least one of the following
two conditions are met.
1. The group containing aj has no other column with a nonzero in row i.
2. The group containing ai has no other column with a nonzero in row j.
Let G(A) = (V, E) be the adjacency graph of A, F od ⊆ E correspond tothe off-diagonal elements in S, and F d correspond to the diagonal elements
in S, i.e., Fd = {(u, u) : u ∈ U } where U ⊆ V . Let F = Fod ∪Fd . A mapping
φ : V → {0, 1, 2, . . . , p} is a distance- 32 coloring of G restricted to F if thefollowing conditions hold.
1. For every (u, u) ∈ Fd ,1.1. φ(u) 6= 0, and1.2. for every (u, v) ∈ E, φ(u) 6= φ(v).
2. For every (v, w) ∈ Fod ,2.1. φ(v) 6= φ(w), and2.2. at least one of the following two conditions holds:
24

2.2.1. φ(v) 6= 0 and for every path (v, w, x), φ(v) 6= φ(x) or
2.2.2. φ(w) 6= 0 and for every path (u, v, w), φ(u) 6= φ(w).
Theorem 4.2 The mapping φ is a distance- 23 p-coloring of G(A) restricted
to F if and only if φ induces a column partition {C 1 , . . . , Cp } symmetricallyconsistent with a direct determination of S.
Proof: Assume that φ is a distance- 23 p-coloring of G(A) restricted to F .
We show that the groups {C1 , . . . , Cp } where Cα = {aj : φ(aj ) = α}, 1 ≤α ≤ p, constitute a partition of a subset of the columns of A symmetricallyconsistent with a direct determination of S.
By coloring conditions 1 and 2, for every a ij ∈ S ((ai , aj ) ∈ F ), at least
one of ai or aj has a nonzero color and hence is involved in the partition{C1 , . . . , Cp }. Let aij ∈ S be a diagonal entry (i = j). Then, coloring
condition 1 ensures that φ(ai ) 6= 0 and that φ(ai ) 6= φ(ak ) for every (ai , ak ) ∈E. Thus, by construction, column ai belongs to group Cφ(ai ) and no column
ak , k 6= i with aik 6= 0 is in Cφ(ai ) . This clearly satisfies the partition
condition. Let aij ∈ S now be an off-diagonal entry (i 6= j). Assume
without loss of generality that φ(a j ) 6= 0. By condition 2.1, φ(aj ) 6= φ(ai ).
By condition 2.2.2, there is no path (a k , ai , aj ) in G(A), for any k 6= i, j suchthat φ(aj ) = φ(ak ). The last two statements together imply that column
aj belongs to group Cφ(aj ) and that no column ak , k 6= j with aik 6= 0 is
in Cφ(aj ) . This satisfies partition condition 1. A similar argument applies
to the case where φ(ai ) 6= 0 which implies the satisfaction of the alternative
partition condition.
To prove the converse, assume that the partition C = {C 1 , . . . , Cp } issymmetrically consistent with a direct determination of S. Construct acoloring φ of G(A) as follows. φ(aj ) = α if aj ∈ Cα , and φ(aj ) = 0 if
aj 6∈ C. We claim that φ is a distance- 23 p-coloring of G(A) restricted to F .
Consider a diagonal element aii ∈ S. The partition conditions ensurethat ai is in some group Cα0 and that there is no column ak ∈ Cα0 , k 6= isuch that aik 6= 0. Thus, by construction, φ(ai ) 6= 0 and φ(ai ) 6= φ(ak )for every (ai , ak ) ∈ E, satisfying coloring condition 1. Consider now thecase where aij ∈ S is an off-diagonal element. First, observe that since all
diagonal elements are nonzero, ai and aj cannot belong to the same group.
Thus φ(ai ) 6= φ(aj ), satisfying coloring condition 2.1. Second, observe that
there are two possibilities by which the partitioning conditions have beensatisfied. We consider only one of these; the second one can be treated ina similar manner. Suppose aj belongs to some group Cα0 and that there isno other column ak ∈ Cα0 , k 6= j such that aik 6= 0. Thus, by construction,25

φ(aj ) 6= 0 and φ(aj ) 6= φ(ak ) for every path (ak , ai , aj ) in G(A), satisfying
coloring condition 2.2.2.
2
A special case of Theorem 4.2 is the problem of estimating only the
diagonal elements of A, i.e., Fd = {(v, v) : v ∈ V } and Fod = ∅. For thisproblem, condition 1 is the only applicable condition, and states that a
distance-1 coloring of G(A) is sufficient.

4.3

Nonsymmetric matrix, Two-dimensional partition

Let A ∈ Rm×n be a nonsymmetric matrix, and S denote the set of nonzero
elements of A required to be estimated. A bipartition (Π C , ΠR ) of a subsetof the columns and rows of A is consistent with a direct determination of S
if for every aij ∈ S at least one of the following conditions are met.
1. The group (in ΠC ) containing column j has no other column having anonzero in row i.
2. The group (in ΠR ) containing row i has no other row having a nonzero
in column j.
Let Gb (A) = (V1 , V2 , E), and F ⊆ E correspond to the elements in S. A
mapping φ : [V1 , V2 ] → {0, 1, . . . , p} is said to be a distance- 23 bicoloring of
Gb restricted to F if the following conditions are met.
1. Vertices in V1 and V2 receive disjoint colors, except for color 0; i.e., forevery u ∈ V1 and v ∈ V2 , either φ(u) 6= φ(v) or φ(u) = φ(v) = 0.
2. At least one endpoint of an edge in F receives a nonzero color; i.e., forevery (v, w) ∈ F , φ(v) 6= 0 or φ(w) 6= 0.
3. For every (v, w) ∈ F ,3.1. if φ(v) = 0, then, for every path (u, v, w), φ(u) 6= φ(w),3.2. if φ(w) = 0, then, for every path (v, w, x), φ(v) 6= φ(x),3.3. if φ(v) 6= 0 and φ(w) 6= 0, then for every path (u, v, w, x), either
φ(u) 6= φ(w) or φ(v) 6= φ(x).
Theorem 4.3 The mapping φ is a distance- 23 p-bicoloring of Gb restricted
to F if and only if φ induces a bipartition (Π C , ΠR ), |ΠC | + |ΠR | = p,consistent with a direct determination of S.
Proof: Let the construction of a partition given a coloring, and vice-versa,
be done in a similar manner as in the proof of Theorem 4.1.
26

Assume that φ is a distance- 32 bicoloring of Gb (A) restricted to F . Let
the induced bipartition be (ΠC , ΠR ). Coloring condition 1 implies that(ΠC , ΠR ) is a bipartition. By condition 2, for every a ij ∈ S, either aj ∈ ΠCor ri ∈ ΠR (or both). Assume now that (ΠC , ΠR ) is not consistent with adirect determination of S. This occurs only if one of the following cases holdfor any aij ∈ S:
• φ(ri ) = 0, φ(aj ) 6= 0 and there exists a column ak , k 6= j with aik 6= 0
such that φ(aj ) = φ(ak ). But this contradicts coloring condition 3.1,
and hence cannot occur.
• φ(aj ) = 0, φ(ri ) 6= 0 and there exists a row rl , l 6= i with alj 6= 0 suchthat φ(ri ) = φ(rl ). But this contradicts coloring condition 3.2, andhence cannot occur.
• φ(ri ) 6= 0, φ(aj ) 6= 0, and there exist column ak , k 6= j with aik 6= 0,
and row rl , l 6= i with alj 6= 0 such that φ(aj ) = φ(ak ) and φ(ri ) =φ(rl ). But this contradicts coloring condition 3.3, and hence cannot
occur.
Hence, (ΠC , ΠR ) is consistent with a direct determination of S.
Conversely, assume that (ΠC , ΠR ) is a bipartition consistent with a direct
determination of S. Clearly, the constructed coloring φ satisfies conditions
1 and 2. To complete the proof, we show that φ also satisfies condition 3.
Assume that φ violates condition 3. Then one of the following cases musthave happened:• There exists a path (ak , ri , aj ) for some (ri , aj ) ∈ F such that φ(ri ) =0 and φ(aj ) = φ(ak ). But this implies that element aij cannot be
determined directly, contradicting the assumption that (Π C , ΠR ) is
consistent with a direct determination of S.
• There exists a path (ri , aj , rl ) for some (ri , aj ) ∈ F such that φ(aj ) = 0
and φ(ri ) = φ(rl ). Again this implies that element a ij cannot be
determined directly, a contradiction of our assumption.
• There exists a path (ak , ri , aj , rl ) for some (ri , aj ) ∈ F such thatφ(ri ) = φ(rl ) 6= 0 and φ(aj ) = φ(ak ) 6= 0. But this implies thatelement aij cannot be determined directly, contradicting the assumption.
2

27

a
a
b
a

b

c

d

e

X
X

h

1

c

X X X

2

d

f

X
X

g

X

3

e

4

f

1

1

2

1

f
2

1

2
2

X X X X

g

1
g

d 1

X X

X X

1

h
2

2

1 b
1

e
1

1

1

c
h

Figure 4: A matrix and its bipartite and intersection graph, resp.

5

Column Intersection vs. Bipartite Graph Formulation

In this section we compare our bipartite graph based formulations withformulations based on a column intersection graph.
Figure 4 depicts a matrix A, the corresponding bipartite graph G b (A) =(V1 , V2 , E), and the column intersection graph G c (A) = (V2 , E 0 ). Note that
we have augmented the graph Gc (A) with edge weights—w(ai , aj ) is the sizeof the intersection of the sets represented by vertices a i and aj . In terms of
matrix A, w(ai , aj ) is the total number of rows where both columns a i andaj have nonzero entries. In the bipartite graph G b (A), it corresponds to thenumber of common neighbors of vertices a i and aj .
In Sections 5.1 to 5.4, we compare and contrast the two graph formulations in terms of flexibility, storage space requirement, ease of graph construction, and use of existing software.

5.1

Flexibility

Notice that the column intersection graph is a ‘compressed’ representationof the structure of the underlying matrix. Clearly, some information is lostin the compression process. In particular, given an edge in G c between twocolumns, we cannot determine the row at which they share nonzero entries.
By contrast, the bipartite graph is an equivalent representation of the structure of the matrix. This provides flexibility. As an illustration, notice thatthe bipartite graph can be used in a one-dimensional (row or column) partition as well as a two-dimensional (combined row and column) partition.
The column intersection graph, on the other hand, is applicable only to a28

column partition. Moreover, the graph formulations for partial matrix estimation problems would have required new types of ‘intersection’ graphs tobe defined. In general, the advantage of the bipartite graph representationis that the representation is decoupled from the eventual technique to be
employed and the matrix entries to be determined.

5.2

Storage Space Requirement

Although Lemma 3.18 correlates the bipartite graph of a matrix with its
column intersection graph, one cannot immediately deduce that one graphis denser than the other. The density of the respective graph namely depends
on the structure of the matrix. Here, we make a rough analysis to show thatfor sparse matrices of practical interest, the column intersection graph islikely to be denser than the bipartite counterpart.
Given a matrix A, the graph Gb (A) = (V1 , V2 , E) and the weighted
graph Gc (A) = (V2 , E 0 ), for a vertex u in Gb , let N1 (u) = {v : (u, v) ∈ E},
d1 (u) = |N1 (u)|, and let the average degree-1 in the sets V 1 and V2 of Gbbe δ 1 (V1 ) and δ 1 (V2 ), respectively. Further, let w denote the average edge
weight in Gc . Then,Xe∈E 0w(e) =1 X X·(d1 (v) − 1)
2 u∈V v∈N (u)211· |V2 | · δ 1 (V2 ) · (δ 1 (V1 ) − 1)
|E 0 | · w ≈2
1=· |E| · (δ 1 (V1 ) − 1)2
δ 1 (V1 ) − 1)|E 0 | = |E| · (2w1 )−1Therefore, as long as δ1 (V2w
> 1, the column intersection graph is likely tohave more edges (and hence requires more storage space) than the bipartite
graph of the matrix.

5.3

Ease of Construction

The sparsity structure of matrix A directly, without any further computation, gives the corresponding bipartite graph G b (A). In principle, the data
structure used to represent A can be used for implementing algorithms thatuse Gb (A). By contrast, Gc (A) has to be computed. As the following Lemma
states, the time required for the computation of G c (A) is proportional to
the number of edges in Gc (A).

29

Lemma 5.1 Given a graph Gb (A) = (V1 , V2 , E), the time required for constructing Gc (A) is Tconst = O(|V2 |(δ 1 (V2 )(δ 1 (V1 ) − 1)).
It should however be noted that once the graph G c (A) is computed, asubsequent distance-1 coloring of G c (A) can be done faster than a distance-2
coloring of Gb (A). As we shall show in Section 6, the overall time requiredfor constructing and distance-1 coloring of G c (A) is of the same order as thetime required for a direct distance-2 coloring of G b (A).

5.4

Use of Existing Software

Serial program packages that implement various practically effective distance1 coloring heuristics exist [5, 6]. For matrix partitioning problems where acolumn intersection graph based formulation can be applied, these packages
can be readily used. On the other hand, since distance-2 coloring is a prototypical model in our context, efficient programs, including parallel ones,
for the distance-2 coloring problem need to be developed.

6

Distance-k Coloring Algorithms

In this section we present several fast algorithms for the various coloring
problems considered in this paper. The algorithm design is greatly simplified
by taking distance-2 coloring as a generic starting point.
For any fixed integer k ≥ 1, DkGCP is NP-hard [21]. A proof sketch
showing that the problem of finding a ( 32 , p)-coloring with the minimum p
is NP-hard, even if the graph is bipartite, is given in [8]. The problem offinding a distance- 23 bicoloring with the fewest number of colors is also NPhard [9]. Further, finding an acyclic p-coloring with the least value of p is
NP-hard [4].
Since all the graph coloring problems of our concern are NP-hard, inpractical applications, we are bound to rely on using approximation algorithms or heuristics. An algorithm A is said to be a γ-approximation algorithm for a minimization problem if its runtime is polynomial in the input
size and if for every problem instance I with an optimal solution OP T (I),A(I)
the solution A(I) output by A is such that OP
T (I) ≤ γ. The approximation
ratio γ ≥ 1, and the goal is to make γ as close to unity as possible. If
no such guarantee can be given for the quality of an approximate solution
obtained by a polynomial time algorithm, the algorithm is usually referred
to as a heuristic.
In the case of distance-1 coloring, there exist several, practically effective
heuristics [7]. In this section we show that some of the ideas used in the30

distance-1 coloring heuristics can be adapted to the cases considered in thispaper by extending the notion of neighborhood. The algorithms we present
are greedy in nature, i.e., the vertices of a graph are processed in some orderand at each step a decision that looks best at the moment (and that will
not be reversed later) is made.
In Section 6.1, we present a generic greedy distance-2 coloring algorithm
and give a detailed analysis of its performance both in terms of computation
time and number of colors used. In Sections 6.2 to 6.5, adaptations of
this algorithm, tailored to the various coloring problems of our concern arepresented.
Notations Some notations used in the rest of this section are first inorder. Recall that for a vertex u in a graph G, a vertex w 6= u is a distancek neighbor of u if the shortest path connecting u and w has length ≤ k.
Let Nk (u) = {w : w is a distance-k neighbor of u}. Let d k (u) = |Nk (u)|denote the degree-k of u; ∆ be the maximum degree-1 in G; and δ k =1 Pu∈V dk (u) denote the average degree-k in G.
|V |Further, in a bipartite graph Gb = (V1 , V2 , E), let the maximum degree-1
in the vertex sets V1 and V2 be denoted by ∆(V1 ) and ∆(V2 ) , respectively.
Similarly, let the average degree-k in the sets V 1 and V2 be denoted by δ k (V1 )and δ k (V2 ), respectively.

6.1

Distance-2 Coloring Algorithms

A simple approach for an approximate distance-2 coloring of a graph G =(V, E) is to visit the vertices in some order, each time assigning a vertex the smallest color that is not used by any of its distance-2 neighbors.
Note that the degree-2 of a vertex u in G is bounded by ∆ 2 , i.e., d2 (u) ≤P
2w∈N1 (u) d1 (w) ≤ ∆ · d1 (u) ≤ ∆ . Thus, since the vertices in G can always
be distance-2 colored trivially using |V | different colors, it is always possible
to color a vertex using a value from the set {1, 2, . . . , min{∆ 2 + 1, |V |}}.
Algorithm GreedyD2Coloring, outlined below, uses this as it colors the vertices of the graph in an arbitrary order. In the algorithm, color(v) is thecolor assigned to vertex v and forbiddenColors is a vector of size C max =min{∆2 + 1, |V |} used to mark the colors that cannot be assigned to aparticular vertex. Specifically, forbiddenColors(c) = v indicates that color ccannot be assigned to vertex v.
Lemma 6.1 GreedyD2Coloring finds a distance-2 coloring in time O(|V |δ 2 ).

31

GreedyD2Coloring(G = (V, E))for each vi ∈ V dofor each colored vertex u ∈ N2 (vi ) do
forbiddenColors(color(u)) = vi
end-forcolor(vi ) = min{c : forbiddenColors(c) 6= vi }end-for

Proof: We first show correctness. In step i of the algorithm, the color usedby each of the distance-2 neighbors of vertex v i is marked (using vi ) in thevector forbiddenColors. Thus, at the end of the inner for loop, the set ofcolors that are allowed for vertex v i is the set of indices in forbiddenColors
where the mark used is different from v i . The minimum value in this set
is thus the smallest allowable color for vertex v i . Notice that the vector
forbiddenColors does not need to be initialized at every step as the marker
vi is used only in step i.
Turning to complexity, note that marking the forbidden colors at step iof the algorithm takes O(d2 (vi )) time. Finding the smallest allowable color
to vi can be done within the same order of time by scanning forbiddenColors
sequentially until the first index c where a value other than v i is stored isP
2found. The total time is thus proportional to v∈V d2 (v) = O(|V |δ 2 ).
We now analyze the quality of the solution provided by GreedyD2Coloring.
Let the number of colors used by GreedyD2Coloring on a graph G = (V, E)
be χ2 greedy (G). Then recalling the lower bound given in Lemma 3.19, we
get the following theorem and its corollary.
Theorem 6.2 ∆ + 1 ≤ χ2 (G) ≤ χ2 greedy (G) ≤ min{∆2 + 1, |V |}.
p

Corollary 6.3 GreedyD2Coloring is an O( |V |)-approximation algorithm.
1
Proof: The approximation ratio γ is at most ∆+1· min{∆2 + 1, |V |}. There
are two possibilities to consider. In the first case ∆ 2 + 1 < |V |. This impliespp
2 +1= O(∆) = O( |V |). In the second case
∆ = O( |V |) and γ = ∆∆+1|V |= O( |V |).
|V | < ∆2 + 1. This implies ∆ = Ω( |V |) and γ = ∆+1

p

p

2

Note that for practical problems, such as problems that arise in solving PDEs using good finite element discretizations, ∆ 2 + 1  |V |, making
GreedyD2Coloring an O(∆)-approximation algorithm.
The actual number of colors used by GreedyD2Coloring depends on theorder in which the vertices are visited. In GreedyD2Coloring, an arbitrary32

ordering is assumed. A solution with fewer number of colors can be expected
if a more elaborate ordering criterion is used. For example, the ideas inlargest degree first and incidence degree ordering for distance-1 coloring [7]can be adapted to the distance-2 coloring case.
As a final remark on the complexity of GreedyD2Coloring, we show thatthe algorithm runs in linear time in the number of vertices for certain sparsePgraphs. Let δ2 = |V1 | u∈V d1 (u)2 and let the standard deviation of degree-1
in G = (V, E) be given by
σ2 =

1 X
2
(d1 (v) − δ 1 ) .
|V | v∈V

Then,|V |σ 2 =X

d1 (v)2 +X

2

δ1 − 2X

d1 (v)δ 1
v∈V
v∈V
X
2
|V |δ2 + |V |δ 1 − 2δ 1
d1 (v)v∈V2
2|V |δ2 + |V |δ 1 − 2|V |δ 1
2
|V |δ2 − |V |δ 1
v∈V

=
=
=
Rewriting we get,2δ2 = δ 1 + σ 2 .
Noting that δ2 ≥ δ 2 , we get the following corollary to Lemma 6.1.
2

Corollary 6.4 GreedyD2Coloring has time complexity O(|V |(δ 1 + σ 2 )).
Since δ 1 = 2|E|
|V | , the complexity expression in Corollary 6.4 reduces to2

O( |E||V | ) for graphs where σ < δ 1 . In particular, for sparse graphs, where|E| = O(|V |), the time complexity of GreedyD2Coloring becomes O(|V |).

6.2

Partial Distance-2 Coloring Algorithms

Here we modify GreedyD2Coloring slightly to make it suitable for solving thepartial distance-2 coloring problem GCP1 (our graph formulation of MPP1).
For any vertex v ∈ V2 , the number of vertices at distance exactly twounits from v is at most ∆(V2 )(∆(V1 ) − 1). Thus, vertex v can always be assigned a color from the set {1, 2, . . . , C max }, where Cmax = min{∆(V2 )(∆(V1 )−1)) + 1, |V2 |}. In Algorithm GreedyPartialD2Coloring, given below, the vector
forbiddenColors is of size Cmax .
The following result is straightforward.

33

GreedyPartialD2Coloring(Gb = (V1 , V2 , E))for each v ∈ V2 dofor each u ∈ N1 (v) dofor each colored vertex w ∈ N1 (u) do
forbiddenColors(color(w)) = v
end-forend-forcolor(v) = min{c : forbiddenColors(c) 6= v}
end-for

Lemma 6.5 GreedyPartialD2Coloring has time complexity
O(|V2 |δ 1 (V2 )(δ 1 (V1 ) − 1)).
As stated earlier, a distance-1 coloring formulation for MPP1 was provided in [7] using the column intersection graph. From Lemmas 5.1 and 6.5
and noting that greedy distance-1 graph coloring is linear in the number ofedges, it follows that the time required for the construction of the column
intersection graph plus the computation of a (greedy) distance-1 coloring is
asymptotically the same as the time required for the direct computation ofa (greedy) partial distance-2 coloring. This means the two formulations are(asymptotically) comparable in terms of overall computation time.

6.3

Distance- 32 Coloring Algorithms

Recall that a distance- 23 coloring, which was used to model MPP2, is adistance-1 coloring where every path of length three uses at least three colors.
We propose two algorithms for this problem. In finding a valid color to assign
a vertex, the first algorithm visits the distance-3 neighbors of the vertex
while the second algorithm visits only the distance-2 neighbors. In bothalgorithms the vector forbiddenColors is of size C max = min{∆2 + 1, |V |}.
GreedyD 23 ColoringAlg1 outlines the first algorithm.
Figure 5 graphically shows the decision made during one of the |V | stepsof GreedyD 23 ColoringAlg1. The root of the tree corresponds to the vertex v
to be colored at the current step. The neighbors of v that are one, two,and three edges away are represented by the nodes at level u, w, and x,
respectively. Each tree node corresponds to many vertices of the input
graph. A darkly shaded node signifies that the vertex is already colored.
The forbidden colors are marked by an f and ‘?’ indicates that whether
the color is forbidden or not depends on the color used at level x. The
correspondence between the figure and Lines 1, 2 and 3 of the algorithm is34

GreedyD 32 ColoringAlg1(G = (V, E))for each v ∈ V dofor each u ∈ N1 (v) do
if u is colored(1)
forbiddenColors(color(u)) = v
for each colored vertex w ∈ N1 (u) do
if u is not colored(2)
forbiddenColors(color(w)) = v
else
for each colored vertex x ∈ N1 (w), x 6= u do
if (color(x) == color (u))(3)
forbiddenColors(color(w)) = v
break
end-ifend-forend-ifend-forend-forcolor(v) = min{c : forbiddenColors(c) 6= v}
end-for

v

f
u
?

f

w

x

Figure 5: Visualizing a step in GreedyD 32 ColoringAlg1

35

obvious.
Notice that in Line 2 of the algorithm, the color of the vertex w in thepath (v, u, w) where u is not yet colored is forbidden for vertex v. Lateron, when vertex u is colored, the test in Line 1 ensures that u gets a color
different from both v and w, making the path use three different colors.
Had the requirement in Line 2 not been imposed, a situation in which apath (v, u, w, x) is two-colored could arise. Thus from its construction, theoutput of GreedyD 32 ColoringAlg1 is a valid distance- 32 coloring. The amountof work done in each step of the algorithm is proportional to d 3 (v). Thus
we get the following result.
Lemma 6.6 GreedyD 23 ColoringAlg1 finds a distance- 23 coloring in timeO(|V |δ 3 ).
We shall now present the second distance- 32 coloring algorithm in which
the coloring at each step is obtained by considering only the distance-2
neighborhood (in contrast to distance-3 neighborhood of the previous case).
The idea behind the algorithm (formulated in terms of matrices) was firstsuggested by Powell and Toint [26].
Recall that distance- 32 coloring is a relaxed distance-2 coloring. As an
illustration, suppose v, u, w, x is a path in a graph. A coloring φ in which
φ(v) = φ(w) = 2, φ(u) = 1 and φ(x) = 3 is a valid distance- 23 (but notdistance-2) coloring on this path.
One way of relaxing the requirement for distance-2 coloring so as toobtain a distance- 23 coloring is to let two vertices at distance of exactly twounits from each other share a color as long as the vertex in between them hasa color of lower value. More precisely, let v, u, w be a path in G and suppose
v and u are colored and we want to determine the color of w. We allow
φ(w) to be equal to φ(v) as long as φ(u) < φ(v). To see that this coloring
can always be extended to yield a valid distance- 32 coloring, consider thepath v, u, w, x, an extension of path v, u, w in one direction. Now, sinceφ(w) = φ(v) > φ(u), we cannot let φ(x) be equal to φ(u). Obviously, φ(x)should be different from φ(w), otherwise it will not be a valid distance-1
coloring. Thus the path v, u, w, x uses three colors, φ is a distance-1 coloring
and therefore it is a valid distance- 23 coloring. The algorithm that makes
use of this idea is given in GreedyD 23 ColoringAlg2.
Clearly, the runtime of GreedyD 32 ColoringAlg2 is O(|V |δ 2 ). Notice, however, that GreedyD 23 ColoringAlg1 may use smaller number of colors than
GreedyD 23 ColoringAlg2. For instance, Figure 6 shows an example where thefirst algorithm uses three colors in coloring the vertices in their alphabeticalorder while the second one uses four in doing the same.
36

GreedyD 32 ColoringAlg2(G = (V, E))for each v ∈ V dofor each u ∈ N1 (v) do
if u is coloredforbiddenColors(color(u)) = v
for each colored vertex w ∈ N1 (u) do
if u is not coloredforbiddenColors(color(w)) = v
else
if (color(w) < color(u))
forbiddenColors(color(w)) = v
end-ifend-forend-forcolor(v) = min{c : forbiddenColors(c) 6= v}
end-for
2

2

b

b

1

3
a

c

1

1

d

3
a

c

4
d

Figure 6: GreedyD 23 ColoringAlg1 vs. Alg2

6.4

Acyclic Coloring Algorithms

Here we present an effective algorithm called GreedyAcyclicColoring to find
an approximate solution for the acyclic coloring problem (GCP3).
The basic idea in our algorithm is to detect and ‘break’ a two-coloredcycle while an otherwise distance-1 coloring of the graph proceeds. Specifically, the algorithm colors the vertices of a graph G = (V, E) while making
a Depth First Search (DFS) traversal. Recall that a back-edge in the DFStree of an undirected graph defines a unique cycle. For more information on
DFS, refer to the book [11].
In GreedyAcyclicColoring, the DFS-tree T (G) of the graph G is implicitly
constructed as the algorithm proceeds. Let φ(v) denote the color of vertex
v and s(v) denote the order in which v is first visited in the DFS traversalof G (1 ≤ s(v) ≤ |V |). The root r of T (G) has s(r) = 1. Further, letp(v) be a pointer to the parent of v in T (G), and l(v) be a pointer to thelowest ancestor of v in T (G) such that φ(l(v)) 6= φ(v) and φ(l(v)) 6= φ(p(v)).
37

The latter pointer will be used in the detection of a two-colored cycle. In
particular, the algorithm proceeds in such a way that the path in T (G) fromv up to, but not including, l(v) is two-colored. At the beginning of thealgorithm, for every vertex u, l(u) is set to point to null.
Consider the step in GreedyAcyclicColoring where vertex v is first visited.
To start with, v is assigned the smallest color different from all of its distance1 neighbors in G, including its parent p(v) in T (G). If there exists a backedge b = (v, w) in the current T (G) such that s(l(p(v))) < s(w) and φ(v) =φ(p(p(v))), then this implies that the cycle corresponding to b is two-colored(see Figure 7 which shows a partial view of the DFS-tree at the step wherevertex v is to be colored). To break the cycle, v is assigned a new color—thesmallest color different from φ(p(v)) and φ(p(p(v)))—and l(v) is set to pointto p(p(v)). Otherwise, if no such back-edge exists, the color of v is declared
final and l(v) is updated in the following manner. If φ(v) 6= φ(p(p(v))), thenl(v) = p(p(v)); otherwise l(v) = l(p(v)).
T(G)
l(p(v))

w?

v

Figure 7: Visualizing a step in GreedyAcyclicColoring
Notice that the amount of work done in each DFS visit of a vertex v in thegraph is proportional to the degree-1 of v. This makes GreedyAcyclicColoringan O(|E|)-time algorithm.
We note that Coleman and Cai [4] have proposed an algorithm for theacyclic coloring problem. The idea in their algorithm is to first transforma given graph G = (V, E) to a ‘completed’ graph G 0 = (V, E 0 ) such that a
distance-1 coloring of G0 is equivalent to an acyclic coloring of G, and thenuse a known distance-1 coloring heuristic on G 0 . The construction of G0 isdone in the following way. Start by setting E 0 = E; visit the vertices in V ina predefined order ; at each step i, if vertex v i is adjacent to vertices vj and38

vk such that both vj and vk are ordered before vi then add the edge (vj , vk )to E 0 .
Our approach differs from that of Coleman and Cai in at least two ways.
First, the graph G0 used in the latter approach may require substantiallymore storage space than the original graph G used in our approach. Second,an edge in E 0 \ E in the latter approach may actually be redundant. For
example, a distance-1 coloring of an odd-length cycle in G uses at least threecolors and hence is a valid acyclic coloring whereas the Coleman and Cai
approach adds one redundant edge to the cycle.

6.5

Bicoloring Algorithms

Here we consider problems GCP4 and GCP5, introduced in Sections 3.4 and
3.5, respectively.
Recall that in a distance- 32 p-bicoloring, some of the vertices are assigned
the neutral color 0. We make the following crucial observation which helps
us identify a possible set of such vertices. The observation is a direct consequence of Condition 2 of Definition 3.8.
Observation 6.7 Let Gb = (V1 , V2 , E) be a bipartite graph andφ : [V1 , V2 ] → {0, 1, . . . , p} be a distance- 23 bicoloring of Gb . Then,• the set C = {v : φ(v) 6= 0} is a vertex cover in G b , and• the set I = {v : φ(v) = 0} is an independent set in G b .
One consequence of Observation 6.7 is that |I| + |C| = |V 1 | + |V2 |. Thus,
minimizing the cardinality of the vertex cover C corresponds to maximizing
the cardinality of the independent set I.
6.5.1

An algorithm for GCP4

In light of Observation 6.7, we suggest GCP4Algorithm as a scheme for solving
the coloring problem GCP4.
GCP4Algorithm(Gb = (V1 , V2 , E))
1. Find a vertex cover C in Gb .
2. Assign the vertices in the set I = (V 1 ∪ V2 ) \ C the color 0.
3. Color the vertices in C such that the result is a distance-3/2 bicoloringof Gb .

39

In Step 1 of GCP4Algorithm, any vertex cover can be used. However,the choice of the vertex cover affects the subsequent coloring in Step 3, bothin terms of number of colors used and coloring time spent. To reduce thecoloring time in Step 3, the size of the vertex cover should be minimized.
Furthermore, minimizing the potential number of colors to be used imposes
an additional requirement: the vertex cover should include those vertices
from V1 and V2 with relatively high number of distance-1 neighbors. (Recall
the introductory discussion in Section 2.2 used to motivate the need fora two-dimensional partition: matrices with a few dense rows and columns
benefit from a two-dimensional partition.)
In a bipartite graph, a minimum cardinality vertex cover can be obtained
via finding a maximum matching in polynomial time [22, 28]. In fact, it canbe computed practically in effectively linear time in the number of edges[24].
Once Steps 1 and 2 are carried out, Step 3 can be done by a suitable
adaptation of GreedyD 32 ColoringAlg1. GreedyD 32 BiColoring, given here only
pictorially, is such an adaptation. One of the differences between the coloring
and bicoloring algorithms is that in the latter case, two disjoint set of colors
are used in coloring the vertices in V 1 and V2 of the bipartite graph Gb =(V1 , V2 , E). Another difference is that at a step of the bicoloring algorithm
where v is colored, a vertex within the distance-3 neighborhood of v may be
one of three types: it is colored with a positive value, it is colored with 0, or
it is not yet colored. The choice of color for v thus needs to consider thesethree options.

v

u
f

f

?
w

x
Figure 8: Visualizing a step in GreedyD 32 BiColoring

40

Figure 8 shows a visual presentation of the ith step of GreedyD 32 BiColoring.
Note the similarity with Figure 5. Notice also that, in choosing a color for
vertex v in Gb = (V1 , V2 , E), since the colors for vertices in V 1 and V2 are
drawn from two disjoint sets, we need only consider colors of vertices twoedges away from v. In the figure, dark shaded nodes correspond to coloredvertices, light shaded nodes show vertices with color zero, and unshaded
nodes correspond to uncolored vertices. Observe that the node with color
0 at level u has only two children. The colors of the vertices in the nodes
marked by an f indicate forbidden colors and whether the color at the nodemarked by ‘?’ is forbidden or not depends on the color used at node x: if
φ(u) = φ(x), φ(w) is forbidden, otherwise, it is not.
The time complexity of GreedyD 32 BiColoring is O((|V1 | + |V2 |)δ 3 ), whichis also the overall time complexity of GCP4Algorithm assuming that step 1
is done using a greedy algorithm that is linear in the number of edges.
Notice that a partial distance-2 coloring of G b on V2 is a just special
case of the scheme GCP4Algorithm. To see this, consider the trivial choice ofvertex cover C = V2 in Step 1. This implies that, in Step 2, the vertices inthe set I = V1 will be colored with zero. By Condition 3 of Definition 3.8,
vertices adjacent to a vertex colored with zero are required to be assigned
different colors. Thus, the result is effectively a partial distance-2 coloringof Gb on V2 .
We note that Hossain and Steihaug [17] and Coleman and Verma [9] haveeach proposed an algorithm for GCP4. These algorithms can be interpreted
in light of GCP4Algorithm. The algorithm of Hossain and Steihaug (HS)implicitly finds a vertex cover while the coloring of the graph proceeds.
Using our terminology, the vertices that remain uncolored at the end of theHS-algorithm form an independent set in the graph and can thus be assigned
the neutral color 0.
The algorithm of Coleman and Verma uses a preprocessing step to identify the rows and columns of the underlying matrix that eventually need
to be colored with positive values. The preprocessing step uses a nonstraightforward matrix-based procedure. It appears that the procedure effectively produces a small sized vertex cover (however, this is not clearly statedin the paper). After the preprocessing step, a certain ‘column intersection’graph, adapted to the distance- 32 bicoloring requirements, is constructed to
finally use known distance-1 coloring heuristics on the resulting graph.

41

6.5.2An algorithm for GCP5

Similarly, by virtue of Observation 6.7, the approach we suggest for solving
the acyclic bicoloring problem is given in GCP5Algorithm.
GCP5Algorithm(Gb = (V1 , V2 , E))
1. Find a vertex cover C in Gb .
2. Assign the vertices in the set I = (V 1 ∪ V2 ) \ C the color 0.
3. Color the vertices in C such that the result is an acyclic bicoloring ofGb .

7

Conclusion

We have studied the efficient estimation of sparse Jacobian and Hessian
matrices using FD and AD techniques. We considered methods that rely
on a one-dimensional as well as a two-dimensional partition to be used inan evaluation based on a direct or a substitution scheme. We introduced
partial matrix estimation problems in distinction from full matrix estimation
problems. In doing so, we developed a unified graph theoretic framework tocope with a variety of complex matrix partitioning problems.
At the basis of our graph problem formulations lies a robust graph representation of the sparsity structure of a matrix: a nonsymmetric matrix isrepresented by its bipartite graph and a symmetric matrix by its adjacency
graph.
We showed that the distance-2 graph coloring problem is a generic modelfor the various one-dimensional matrix partitioning problems.
Our unified graph theoretic approach enabled us to provide some fresh
insight into the matrix problems and as a result we developed several simple and effective algorithms. Our emphasis has been on greedy algorithms.
Other algorithmic techniques need to be explored in the future. For example, it could be interesting to find a distance-2 coloring algorithm thatuses asymptotically the same time as the greedy algorithm discussed in thispaper and balances the number of vertices in each color class. Finding arandom color, rather than the smallest color, from an allowable set could be
an idea to consider in this regard.
In the case of two-dimensional partition problems, based on the known
relationship to graph bicoloring, we argued that finding a ‘small’-size vertex
cover as a preprocessing step contributes to making the overall computation
42

more efficient. A more precise characterization of the ‘optimum’ vertex coverrequired is a worthwhile issue.
We have not developed any special algorithms for the restricted coloring
problems arising in partial matrix estimation. The ideas used in our algorithms for the coloring problems in full matrix estimation can be adapted
to the restricted cases by observing the particular coloring conditions.
In general, most of the algorithms in the literature for solving the coloring
problems considered in this paper rely on first transforming the input graph
G = (V, E) to some denser graph G0 = (V, E 0 ), E 0 ⊇ E, such that a distance1 coloring of G0 is equivalent to the particular coloring problem on G. In
contrast, the algorithms proposed in this paper solve the particular coloring
problem directly on G. As has been argued, the main advantages offeredby our approach are the possibility to mix-and-match methods, less storage
space requirement, and ease of developing flexible software.
One of the motivations for the current study has been the need for thedevelopment of parallel algorithms for solving partitioning problems in largescale PDE-constrained optimization contexts. In a recent work [13], we haveshown some parallel algorithms (using the shared-memory programming
model) for the distance-2 and distance- 32 coloring problems. Our results,
theoretical as well as experimental, were promising. We believe that thisstudy lays a foundation for further work on the development and implementation of not only shared-memory but also distributed-memory parallel
algorithms.
Acknowledgement

We thank Trond Steihaug for interesting discussions.

References[1] G. Agnarsson, R. Greenlaw, and M. M. Halldórsson. On powers ofchordal graphs and their colorings. Congress Numerantium, 100:41–65,
2000.
[2] G. Agnarsson and M. M. Halldórsson. Coloring powers of planar graphs.
In 11th. Ann. ACM-SIAM symp. on Discrete Algorithms, pages 654–
662, 2000.
[3] O. Axelsson and U. Nävert. On a graphical package for nonlinear partial differential equation problems. In B. Gilchrist, editor, Proceedingsof IFIP Congress 77, Information Processing, pages 103–108. NorthHolland, 1977.

43

[4] T. F. Coleman and J. Cai. The cyclic coloring problem and estimationof sparse Hessian matrices. SIAM J. Alg. Disc. Meth., 7(2):221–235,
April 1986.
[5] T. F. Coleman, B. Garbow, and J. J. Moré. Software for estimating
sparse Jacobian matrices. ACM Trans. Mathematical Software, 10:329–
347, 1984.
[6] T. F. Coleman, B. Garbow, and J. J. Moré. Software for estimating
sparse Hessian matrices. ACM Trans. Mathematical Software, 11:363–
377, 1985.
[7] T. F. Coleman and J. J. Moré. Estimation of sparse Jacobian matrices
and graph coloring problems. SIAM J. Numer. Anal., 20(1):187–209,
February 1983.
[8] T. F. Coleman and J. J. Moré. Estimation of sparse Hessian matrices
and graph coloring problems. Math. Program., 28:243–270, 1984.
[9] T. F. Coleman and A. Verma. The efficient computation of sparse Jacobian matrices using automatic differentiation. SIAM J. Sci. Comput.,
19(4):1210–1233, July 1998.
[10] G. Corliss, C. Faure, A. Griewank, L. Hascoét, and U. Naumann. Automatic Differentiation of Algorithms: From Simulation to Optimization.
Springer, New York, 2002.
[11] T. H. Cormen, C. E. Leiserson, and R. L. Rivest. Introduction to Algorithms. Massachusetts Institute of Technology, 1990.
[12] A. R. Curtis, M. J. D. Powell, and J. K. Reid. On the estimation ofsparse Jacobian matrices. J. Inst. Math. Appl., 13:117–119, 1974.
[13] A.H. Gebremedhin, F. Manne, and A. Pothen. Parallel distance-kcoloring algorithms for numerical optimization. In B. Monien and
R. Feldmann, editors, Euro-Par 2002 Parallel Processing, volume 2400
of LNCS, pages 912–921. Springer Verlag, 2002.
[14] A. Griewank and G. F. Corliss. Automatic Differentiation of Algorithms: Theory, Implementation, and Application. SIAM, Philadelphia,
1991.
[15] B. Grünbaum. Acyclic colorings of planar graphs. Israel J. Math.,
14:390–408, 1973.
44

[16] S. Hossain. On the computation of sparse Jacobian matrices and Newton steps. Technical Report 146, Department of Informatics, Universityof Bergen, Norway, March 1998.
[17] S. Hossain and T. Steihaug. Computing a sparse Jacobian matrix byrows and columns. Optimization Methods and Software, 10:33–48, 1998.
[18] S. Hossain and T. Steihaug. Reducing the number of AD passes for computing a sparse Jacobian matrix. In G. Corliss, C. Faure, A. Griewank,
L. Hascoét, and U. Naumann, editors, Authomatic Differentiation of Algorithms: From Simulation to Optimization, pages 263 – 270. Springer,
2002.
[19] D. A. Knoll and D. E. Keyes.
Jacobian-free Newton-Krylov
methods: A survey of appoaches and applications. Available at
http://www.math.odu.edu/˜ keyes/.
[20] S. O. Krumke, M. V. Marathe, and S.S. Ravi. Approximation algorithms for channel assignment in radio networks. In Dial M for Mobility,
Dallas, Texas, October 1998.
[21] Y. Lin and S. S. Skiena. Algorithms for square roots of graphs. SIAM
J. Disc. Math., 8:99–118, 1995.
[22] L. Lovász and M.D. Plummer. Matching Theory. North-Holland, Amsterdam, 1986.
[23] S. T. McCormick. Optimal approximation of sparse Hessians and its
equivalence to a graph coloring problem. Math. Program., 26:153–171,
1983.
[24] R. Motwani. Average-case analysis of algorithms for matchings andrelated problems. Journal of the ACM, 41:1329–1356, 1994.
[25] G. N. Newsam and J. D. Ramsdell. Estimation of sparse Jacobianmatrices. SIAM J. Alg. Disc. Meth., 4:404–418, 1983.
[26] M. J. D. Powell and P. L. Toint. On the estimation of sparse Hessian
matrices. SIAM J. Numer. Anal., 16(6):1060–1074, December 1979.
[27] V. V. Vazirani. Approximation Algorithms. Springer, 2001.
[28] L. A. Wolsey. Integer Programming. Wiley-Interscience Publication,
New York, USA, 1998.

45

Parallel Distance-k Coloring Algorithms for
Numerical Optimization
Assefaw Hadish Gebremedhin∗

Fredrik Manne

Alex Pothen‡†

Abstract
Matrix partitioning problems that arise in the efficient estimation ofsparse Jacobians and Hessians can be modeled using variants of graphcoloring problems. In a previous work [7], we argue that distance2 and distance- 23 graph coloring are robust and flexible formulationsof the respective matrix estimation problems. The problem size inlarge-scale optimization contexts makes the matrix estimation phasean expensive part of the entire computation both in terms of execution
time and memory space. Hence, there is a need for both shared- anddistributed-memory parallel algorithms for the stated graph coloring
problems. In the current work, we present the first practical sharedaddress space parallel algorithms for these problems. The main idea in
our algorithms is to randomly partition the vertex set equally among
the available processors, let each processor speculatively color its vertices using information about already colored vertices, detect eventual
conflicts in parallel, and finally re-color conflicting vertices sequentially.
Randomization is also used in the coloring phases to further reduce
conflicts. Our PRAM-analysis shows that the algorithms should givealmost linear speedup for sparse graphs that are large relative to thenumber of processors. Experimental results from our OpenMP implementations on a Cray Origin2000 using various large graphs show thatthe algorithms indeed yield reasonable speedup for modest numbers ofprocessors.

∗

Department of Informatics, University of Bergen, N-5020 Bergen, Norway. {assefaw,
fredrikm}@ii.uib.no†
Computer Science Department, Old Dominion University, Norfolk, VA 23529 USA;
CSRI, Sandia National Labs, Albuquerque NM 87185 USA; and ICASE, NASA LangleyResearch Center, Hampton, VA 23681-2199 USA. pothen@cs.odu.edu‡
This author’s research was supported by NSF grant DMS-9807172, DOE ASCI level-2
subcontract B347882 from Lawrence Livermore National Lab; and by DOE SCIDAC grant
DE-FC02-01ER25476.

1

1Introduction

Numerical optimization algorithms that rely on derivative information often
need to compute the Jacobian or Hessian matrix. Since this is an expensive part of the computation, efficient methods for estimating these matrices
via finite differences (FD) or automatic differentiation (AD) are needed. It
is known that the problem of minimizing the number of function evaluations (or AD passes) required in the computation of these matrices can beformulated as variants of graph coloring problems [1, 2, 3, 9, 13]. The particular coloring problem differs with the optimization context: whether theJacobian or the Hessian matrix is to be computed; whether a direct or asubstitution method is employed; and whether only columns, or only rows,
or both columns and rows are to be used to evaluate the matrix elements. In
addition, the type of coloring problem depends on the kind of graph used torepresent the underlying matrix. In [7], we provide an integrated review ofprevious works in this area and identify the distance-2 (D2) graph coloring
problem as a unifying, generic, and robust formulation. The D2-coloring
problem has also noteworthy applications in other fields such as channel
assignment [11] and facility location problems [14].
Large-scale PDE-constrained optimization problems can be solved onlywith the memory and time resources available on parallel computers. In
these problems, the variables defined on a computational mesh are already
distributed on the processors, and hence parallel coloring algorithms areneeded for computing, for instance, the Jacobian. It turns out that the
problems of efficiently computing the Jacobian and Hessian can be formulated as the D2- and D 32 -coloring problems, respectively. The latter coloring
problem is a relaxed variant of the former, and will be described in Section 2.2. In these formulations, the bipartite graph associated with the rows
and columns of the matrix is used for the Jacobian; in the case of the Hessian matrix, the adjacency graph corresponding to the symmetric matrix isused.
In this paper, we present several new deterministic as well as probabilistic
parallel algorithms for the D2- and D 32 -coloring problems. Our algorithms
are practical and effective, well suited for shared address space programming,and have been implemented in C using OpenMP primitives. We reportresults from experiments conducted on a Cray Origin 2000 using large graphsthat arise in finite element methods and in eigenvalue computations.
In the sequel, we introduce the graph problems in Section 2, present thealgorithms in Section 3, discuss our experimental results in Section 4, andconclude the paper in Section 5.
2

2

Background

2.1

Matrix Partition Problems

An essential component of the efficient estimation of a sparse Jacobian or
Hessian using FD or AD is the problem of finding a suitable partition of
the columns and/or the rows of the matrix. The particular partition chosen
defines a system of equations from which the matrix entries are determined.
A method that utilizes a diagonal system is called a direct method, andone that uses a triangular system is called a substitution method. A direct
method is more restrictive but the computation of matrix entries is straightforward and numerically stable. A substitution method, on the other hand,is less restrictive but it may be subject to approximation difficulties andnumerical instability. Moreover, direct methods offer more parallelism thansubstitution methods. In this paper, we focus on direct methods that use
column partitioning.
A partition of the columns of a nonsymmetric matrix A is said to beconsistent with the direct determination of A if whenever a ij is a non-zeroelement of A then the group containing column j has no other column witha non-zero in row i [1]. Similarly, a partition of the columns of a symmetricmatrix A is called symmetrically consistent with the direct determinationof A if whenever aij is a non-zero element of A then either (i) the groupcontaining column j has no other column with a non-zero in row i, or (ii)the group containing column i has no other column with a non-zero in row
j [2]. From a given (symmetrically) consistent partition {C 1 , C2 , . . . , Cρ } ofthe columns of A, the nonzero entries can be determined with ρ function
evaluations (matrix-vector products).
Thus we have the following two problems of interest. Given the sparsity
structure of a nonsymmetric m × n matrix A, find a consistent partition of
the columns of A with the fewest number of groups. We refer to this problem
as NONSYMCOLPART. The second problem of our interest, which we call
SYMCOLPART, states: given the sparsity structure of a symmetric n × nmatrix A, find a symmetrically consistent partition of the columns of A withthe fewest number of groups.

2.2Graph Problems

In a graph, two distinct vertices are said to be distance-k neighbors if theshortest path connecting them consists of at most k edges. The numberof distance-k neighbors of a vertex u is referred to as the degree-k of u
and is denoted by dk (u). A distance-k ρ-coloring (or (k, ρ)-coloring for3

short) of a graph G = (V, E) is a mapping φ : V → {1, 2, . . . , ρ} suchthat φ(u) 6= φ(v) whenever u and v are distance-k neighbors. We call amapping φ : V → {1, 2, . . . , ρ} a ( 32 , ρ)-coloring of a graph G = (V, E) if φ
is a (1, ρ)-coloring of G and every path containing three edges uses at least
three colors. Notice that a ( 23 , ρ)-coloring is a restricted (1, ρ)-coloring, anda relaxed (2, ρ)-coloring, and hence the name. For instance, consider a pathu, v, w, x in a graph. The assignment 2, 1, 2, 3 to the respective vertices isa valid D1- and D 23 -coloring, but not a valid D2-coloring. The distance-k
graph coloring problem asks for a (k, ρ)-coloring of a graph with the least
possible value of ρ.
Let A be an m × n rectangular matrix with rows r 1 , r2 , . . . , rm andcolumns a1 , a2 , . . . , an . We define the bipartite graph of A as G b (A) =(V1 , V2 , E) where V1 = {r1 , r2 , . . . , rm }, V2 = {a1 , a2 , . . . , an }, and (ri , aj ) ∈E whenever aij 6= 0, for 1 ≤ i ≤ m, 1 ≤ j ≤ n. When A is an n × n symmetric matrix with non-zero diagonal elements, the adjacency graph of A isdefined to be Ga (A) = (V, E) where V = {a1 , a2 , . . . , an } and (ai , aj ) ∈ E
whenever aij , i 6= j, is a non-zero element of A. Note that the non-zerodiagonal elements of A are not explicitly represented by edges in G a (A).
In our work, we rely on the bipartite and adjacency graph representations.
However, in the literature, a nonsymmetric matrix A is often represented byits column intersection graph Gc (A). In this representation, the columns ofA constitute the vertex set, and an edge (a i , aj ) exists whenever the columns
ai and aj have non-zero entries at the same row position (i.e., a i and aj are
not structurally orthogonal). As argued in [7], the bipartite graph representation is more flexible and robust than the ‘compressed’ column intersection
graph.
Coleman and Moré [1] showed that problem NONSYMCOLPART isequivalent to the D1-coloring problem when the matrix is represented byits column intersection graph. We have shown [7] that the same problem
is equivalent to a partial D2-coloring when a bipartite graph is used anddiscussed the relative merits of the two approaches. The word ‘partial’ reflects the fact that only the vertices corresponding to the columns need tobe colored.
McCormick [13] showed that the approximation of a Hessian using a
direct method is equivalent to a D2-coloring on the adjacency graph of thematrix. One drawback of McCormick’s formulation is that it does not exploit
symmetry. Later Coleman and Moré [2] addressed this issue and showed thatthe resulting problem (SYMCOLPART) is equivalent to the D 23 -coloring
problem.

4

3Parallel Coloring Algorithms

The distance-k coloring problem is NP-hard for any fixed integer k ≥ 1[12]. A proof-sketch showing that D 32 -coloring is NP-hard is given in [2].
Furthermore, Lexicographically First ∆ + 1 Coloring (LFC), the polynomial
variant of D1-coloring in which the vertices are given in a predeterminedorder and the question at each step is to assign the vertex the smallest color
not used by any of its neighbors, is P-complete [8]. The practical implicationof this is that designing efficient fine-grained parallel algorithm for LFC ishard.
In practice, greedy sequential D1-coloring heuristics are found to be quite
effective [1]. In a recent work [6], we have shown effective methods of parallelizing such greedy algorithms in a coarse-grained setting. Here, we extend
this work to develop parallel algorithms for the D2- and D 23 -coloring problems.
Jones and Plassmann [10] describe a parallel distributed memory D1coloring algorithm that uses randomization to assign priorities to the vertices, and then colors the vertices in the order determined by the priorities.
There is no speculative coloring in their algorithm. It was reported that thealgorithm slows down as the number of processors is increased. Finocchi etal. [5] suggest a parallel D1-coloring algorithm organized in several rounds;
in each round, currently uncolored vertices are assigned a tentative pseudocolor without consulting their neighbors mapped to other processors; in a
conflict resolution step, a maximal independent set of vertices in each color
class is assigned these colors as final; the remainder of the vertices are uncolored, and the algorithm moves into the next round. However, they do
not give any implementation and we believe that this algorithm incurs too
many rounds, each with its synchronization and communication steps, for
it to be practical on large graphs.
Our algorithm (in its generic form) may be viewed as a compromise
between these algorithms, where we permit speculative coloring, but limit
the number of synchronization steps to two in the whole algorithm. However,
our current algorithms rely on the shared address space programming model;
we will adapt our algorithms to distributed memory programming models
in future work. We are unaware of any previous work on parallel algorithms
for D2- and D 23 -coloring problems.

5

3.1

Generic Greedy Parallel Coloring Algorithm (GGPCA)

The steps of our generic parallel coloring algorithm can be summarized asshown below; refer to [6] for a detailed discussion of the D1-coloring case.
Let G = (V, E) be the input graph and p be the number of processors.
Phase 0 : Partition
Randomly partition V into p equal blocks V 1 . . . Vp . Processor Pi is
responsible for coloring the vertices in block V i .
Phase 1 : Pseudo-color
for i = 1 to p do in parallelfor each u ∈ Vi do
assign the smallest available color to u, paying attention
to already colored vertices (both local and non-local).
—barrier synchronize—Phase 2 : Detect conflicts
for i = 1 to p do in parallelfor each u ∈ Vi docheck whether the color of u is valid. If the colors of u and v
are the same for some ‘neighbor’ v of u thenLi = Li ∪ min{u, v}—barrier synchronize—Phase 3 : Resolve conflicts
Color the vertices in the conflict list L = ∪L i sequentially.
In Phase 0, the vertices are randomly partitioned into p equal blocks eachof which is assigned to some processor. In Phase 1, the processors color thevertices in their respective blocks in parallel. When two ‘neighbor’ vertices
reside on different processors, the two processors could color both simultaneously, possibly assign them the same value, and cause a conflict. The
purpose of Phase 2 is to detect and store any such conflict vertices which
are subsequently re-colored sequentially in Phase 3.

3.2

Simple Distance-2 (SD2) Coloring Algorithm

The meaning of ‘available’ color in GGPCA depends on the required coloring. In the case of a D2-coloring, a vertex is assigned the smallest color
not used by any of its distance-2 neighbors. Let ∆ denote the maximum
degree-1 in the graph G = (V, E). It can easily be verified that, in a D2coloring, a vertex can always be assigned one of the colors from the set{1, 2, . . . , ∆2 + 1}. Moreover, since the distance-1 neighbors of a vertex are6

distance-2 neighbors with each other, the 2-chromatic number, i.e., the least
number of colors required in a D2-coloring, is at least ∆+1. Thus, the greedy
approach is an O(∆)-approximation algorithm. We refer to the variant of
GGPCA that applies to D2-coloring as Algorithm SD2.
Note that the sequential time complexity of greedy D2-coloring is O(∆ 2 |V |).
The following results show that the number of conflicts discovered in Phase
2 of Algorithm SD2 is often small for sparse graphs,
r making the algorithm
scalable when the number of processors p = O(

|V |2∆|E| ).

The proofs are

straightforward extensions of our proofs for the D1-coloring case given in [6]and hence are omitted here. One only needs to observe that the degree-2
and degree-1 of a vertex u are related by d 2 (u) ≤ ∆d1 (u). Let δ = 2|E|/|V |denote the average degree-1 in G.
Lemma 1 The expected number of conflicts created at the end of Phase 1.
of Algorithm SD2 is at most ≈ ∆δ(p−1)
2Theorem 2 On a CREW PRAM, Algorithm SD2 distance-2 colors the input graph consistently in expected time O(∆ 2 ( |Vp | + ∆δp)) using at most
∆2 + 1 colors.
Corollary 3 When p = O(

r

|V |2
∆|E| ), the expected runtime of Algorithm SD2

2

is O( ∆ p|V | ).
The number of conflicts predicted by Lemma 1 is an overestimate. The
analysis assumes that whenever two distance-2 neighbor vertices are colored simultaneously, they are assigned the same color, thereby resulting in a
conflict. However, a more involved probabilistic analysis that takes the distribution of colors used into account may provide a tighter bound. Besides,
the actual number of conflicts in an implementation could be significantly
reduced by choosing a random color from the allowable set, instead of thesmallest one as given in Phase 1 of GGPCA.

3.3

Improved Distance-2 (ID2) Coloring Algorithm

The number of colors used in Algorithm SD2 can be reduced using a ‘tworound-coloring’ strategy. The underlying idea in the D1-coloring case is due
to Culberson [4] and was used in our parallel algorithms for D1-coloring[6]. In Lemma 4 we extend the result to the D2-coloring case; the proof
is basically a reproduction of Culberson’s proof for the distance-1 coloringcase. The greedy sequential algorithm referred to in the lemma is one that7

visits the vertices of a graph G = (V, E) in some order (a permutation of{1, 2, . . . , |V |}), each time assigning a vertex the smallest allowed color. In
particular, the first vertex to be visited is assigned color 1.
Lemma 4 Let φ be a distance-2 coloring of a graph G using α colors, andπ a permutation of the vertices such that if φ(v π(i) ) = φ(vπ(l) ) = c, thenφ(vπ(j) ) = c for i < j < l. Applying the greedy sequential distance-2 coloringalgorithm to G where the vertices have been ordered by π will produce acoloring φ0 using α or fewer colors.
Proof: The proof is a simple induction showing that the first i color classes 1
listed in the permutation will be colored with i or fewer colors. Clearly, thefirst color class listed will be colored with color 1. Suppose some element of
the ith class requires color i + 1. This means that it must have a distance-2
neighbor of color i. But by induction the vertices in the 1st to the (i − 1)thclasses used no more than i − 1 colors. Thus, the conflict has to be with amember of its own color class, but this contradicts the assumption that φ isa valid distance-2 coloring.
2
The idea in Lemma 4 is that if the greedy coloring algorithm is re-applied
on a graph, with the vertices belonging to the same color class (in the originalcoloring) listed consecutively, then the new coloring obtained is better or atleast as good as the original. One ordering (among many) that satisfies thiscondition, with a good potential for reducing the number of colors used, is tolist the vertices consecutively for each color class in the reverse order of theintroduction of the color classes. Based on Lemma 4, we modify Algorithm
SD2 and introduce an additional parallel coloring phase between Phases 1
and 2. Algorithm ID2 below outlines the resulting 4-phase algorithm.
Phases 0 and 1. Same as Ph. 0 and 1 of GGPCA.
(Let s be the number of colors.)
Phase 2. for k = s downto 1 do
Partition ColorClass(k) into p equal blocks V 10 , . . . , Vp0
for i = 1 to p do in parallelfor each u ∈ Vi0 do
assign the smallest available color to vertex u.
—barrier synchronize—Phases 3 and 4. Same as Phases 2 and 3 of GGPCA, respectively.
1

Vertices of the same color constitute a color class.

8

In Phase 2, most of the vertices in a color class are at a distance greaterthan two edges from each other, and the exceptions arise from the conflictvertices colored incorrectly in Phase 1. Since the number of such conflictvertices from Phase 1 is low, the number of conflict vertices at the end of there-coloring phase will be even lower. Phases 3 and 4 are included to detect
and resolve any eventual conflicts not resolved in Phase 2.

3.4

Simple Distance- 32 (SD 23 ) Coloring Algorithm

Recall that a D 32 -coloring is a relaxed D2-coloring (see the example in Section 2.2). One way of relaxing the requirement for D2-coloring in GGPCA
so as to obtain a valid D 32 -coloring is to let two vertices at a distance ofexactly two edges from each other share a color as long as the vertex in
between them is already colored with a (color of) lower value. We refer tothe variant of GGPCA that employs this technique to achieve a distance-3/2
coloring as Algorithm SD 32 .
Note that both Algorithms ID2 and SD 32 take asymptotically the sametime as Algorithm SD2.

3.5

Randomization

The potential scalability of GGPCA depends on the number of conflicts discovered in Phase 2, since these are resolved sequentially in Phase 3. For
dense graphs the number of conflicts could be large enough to destroy thescalability of the algorithm. To overcome this problem, we use randomization as a means for reducing the number of conflicts. In the randomized
variants of Algorithms SD2, SD 32 , and ID2, a vertex is assigned the next
available color with probability q, where 0 < q ≤ 1. The first attempt ismade with the smallest available color, i.e., the color is chosen with probability q. An attempt is said to be successful if the vertex is assigned a color.
If an attempt is not successful, then the next available color is tried withprobability q, and so on, until the vertex gets a color. Algorithms SD2, SD 23 ,and ID2, can be seen as the deterministic variants where q = 1. We refer tothe randomized versions of the respective algorithms as RSD2, RSD 32 , andRID2.
Let u and v be two vertices with the same (infinite) set of allowable colors.
If u and v are colored concurrently, it can be shown that the probability thatu and v get the same color is q/(2 − q). This shows how randomization leads
to a reduction in the number of conflicts; the lower the value of q, the lower
chance for a conflict to arise. It should however be noted that a ‘low’ value

9

Set
Set I (FE)

Problem
mrng2
fe144
m14b

|V |
1,017,253
144,649
214,765|E|
2,015,714
1,074,393
1,679,018

∆
4
26
40

δ
2
4
4

δ
4
15
16

Set II (EV)

ev01
ev02

10,134
19,845

1,318,579
3,353,890

634
749

93
78

260
338

Table 1: Test Graphs: the last three columns list the max., min., and average
degree-1, respectively.
of q may result in an increase in the number of colors used. Choosing theright value for q thus becomes a design issue.

4

Experimental Results and Discussion

Our test bed consists of graphs that arise from finite element methods andfrom eigenvalue computations [6]. Table 1 provides the test graphs’ structural information: columns |V | and |E| give the number of vertices andedges of each graph, respectively, and the maximum, minimum, and average
degree-1 of each graph is given under columns ∆, δ, and δ, respectively.
Table 2 through 7 provide coloring and timing information of the different
algorithms. The number of blocks (processors) is given in column p. Column
χi gives the number of colors used at the end of Phase i of the corresponding
algorithm. The number of conflicts that arise in Phase 1 (or 2) are listed
under the column labeled K (or K2 ). The time in milliseconds requiredby the different phases are listed under T 1 , T2 , T3 , T4 ; column Ttot gives
the total time used. The last two columns display speedup with respect
to two different references: Sseq lists the speedup obtained in comparison
with the runtime of a pure sequential version, i.e., a version with no conflict
detection phase, and Spar displays speedup obtained by taking the runtime
of a parallel algorithm on one processor as a reference.
Since the deterministic algorithms gave acceptable results on the relatively sparse graphs of Set I, the randomized variants were run only on thegraphs from Set II.
The experimental results show that Algorithm SD2 uses many fewer
colors than the bound ∆2 + 1 and that Algorithm ID2 reduces the numberof colors by up to 10% compared to SD2. The advantage of exploiting
symmetry in Problem SYMCOLPART can be seen by comparing the numberof colors used in SD2 and SD 32 ; the latter can be as much as 37% fewer
than the former. In general, the number of colors used in our deterministic
algorithms increases only slightly with increasing p. For the randomized

10

Fig. 1a Simple D2

Fig. 1b Improved D2

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0
p= 1 2 4 6 8 12

1 2 4 6 8 12

Fig. 1c Simple D3/2

Ph. 1
Ph. 2
Ph. 3
Ph. 4

0
p= 1 2 4 6 8 12

Fig. 1d Randomized Algorithms

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0
p= 1 2 4 6 8 12

1 2 4 6 8 120
p=1 2 4 6 8 12 1 2 4 6 8 12 1 2 4 6 8 12

1 2 4 6 8 12

RSD2

RID2

RSD3/2

Figure 1: a-c: Relative performance of deterministic algorithms on graphs
fe144 (left-half on each sub-figure) and ev02 (right-half). Fig. 1d: Relative
performance of randomized algorithms on ev02.
algorithms, the number of colors in some cases even decreases as p increases,
but we think this is a random phenomenon.
Based on the results displayed in Tables 2 through 7, Fig. 1 shows how
the different phases of the algorithms scale as more processors are employed
for two representative graphs: fe144 from Set I and ev02 from Set II. Forgraph fe144, the time used in resolving conflicts sequentially is negligible
compared to the overall time. Moreover, it can be seen that Phases 1 and 2of Algorithms SD2 and SD 32 scale rather well on this graph as the numberof processors is increased. However, Phase 2 of Algorithm ID2 does notscale as well. This is due to the existence of many color classes with few
vertices which entails extra synchronization and communication overhead in
the parallel re-coloring phase.
The picture for the more dense graph ev02 is different: the time elapsedin Phase 3 of Algorithm SD2 is significant and increases as the number ofprocessors is increased. The situation is somewhat better for SD3/2 andID2. Note that ev02 is about 100 times denser than fe144 (where density= |V|E||2 ), and the results in Tables 2 to 7 and Fig. 1 agree well with the results
11

in Corollary 3: when the average degree is high, we lose scalability.
Fig. 1d shows how using probabilistic algorithms solves the problem ofhigh number of conflicts for graph ev02. The improvement in scalability
comes at the expense of increased number of colors used (see Table 2 to 7).
We have experimented using different values for q and found good results
when q = 1/20 for RSD2, and q = 1/6 for RID2 and RSD 32 .
It should be noted that the speedups observed in Fig. 1 are all relative
to the respective parallel algorithm run with p = 1 (see S par in the various
tables). A comparison against a sequential version with no conflict detecting
and resolving phase would yield less speedup (see S seq in the tables). In
particular, relative to a sequential version, the ideal speedup obtained by
Algorithms SD2 and ID2 is roughly 12 p and 23 p, respectively.
Our algorithms in general did not scale well beyond around 16 processors. We believe this is due to, among other things, the relatively high costassociated with non-local physical memory accesses. It would be interesting
to see how this affects the behavior of the algorithms on different parallel
platforms.

5

Conclusion

We have presented several simple and effective parallel approximation algorithms as well as results from OpenMP-implementations for the D2- andD 32 -coloring problems. The number of colors produced by the algorithms in
the case where p = 1 is generally good as it is typically off from the lower
bound ∆ + 1 of SD2 by a factor much less than the approximation ratio ∆.
As more processors are employed, the algorithms provide reasonable speedup
while maintaining the quality of the solution. In general, our deterministic
algorithms seem to be suitable for sparse graphs and the probabilistic variants for more dense graphs. We believe the functionality provided by ouralgorithms is useful for many large-scale optimization codes, where parallel speedups while desirable, are not paramount, as long as running times
for coloring are low relative to the other steps in the optimization computations. The three sources of parallelism in our algorithms – partitioning,
speculation, and randomization – can be exploited in developing distributed
parallel algorithms, but the algorithms would most likely differ significantly
from the shared memory variants presented here.

12References[1] T. F. Coleman and J. J. Moré. Estimation of sparse Jacobian matrices
and graph coloring problems. SIAM J. Numer. Anal., 20(1):187–209,
February 1983.
[2] T. F. Coleman and J. J. Moré. Estimation of sparse Hessian matrices
and graph coloring problems. Math. Program., 28:243–270, 1984.
[3] T. F. Coleman and A. Verma. The efficient computation of sparse Jacobian matrices using automatic differentiation. SIAM J. Sci. Comput.,
19(4):1210–1233, July 1998.
[4] J. C. Culberson. Iterated greedy graph coloring and the difficulty landscape. Technical Report TR 92-07, Department of Computing Science,
The University of Alberta, Edmonton, Alberta, Canada, June 1992.
[5] I. Finocchi, A. Panconesi, and R. Silvestri. Experimental analysis
of simple, distributed vertex coloring algorithms. In Proceedings of
the Thirteenth ACM-SIAM Symposium on Discrete Algorithms (SODA02), San Francisco, CA, 2002.
[6] A. H. Gebremedhin and F. Manne. Scalable parallel graph coloringalgorithms. Concurrency: Pract. Exper., 12:1131–1146, 2000.
[7] A. H. Gebremedhin, F. Manne, and A. Pothen.
Graph coloring in optimization revisited. Technical Report 226, University ofBergen, Dept. of Informatics, Norway, January 2002. Available at:
http://www.ii.uib.no/publikasjoner/texrap/.
[8] R. Greenlaw, H.J. Hoover, and W. L. Ruzzo. Limits to Parallel Computation: P-Completeness Theory. Oxford University Press, New York,1995.
[9] A.K.M S. Hossain and T. Steihaug. Computing a sparse Jacobian matrix by rows and columns. Optimization Methods and Software, 10:33–
48, 1998.
[10] M. T. Jones and P. E. Plassmann. A parallel graph coloring heuristic.
SIAM J. Sci. Comput., 14(3):654–669, May 1993.
[11] S. O. Krumke, M. V. Marathe, and S. S. Ravi. Approximation algorithms for channel assignment in radio networks. In Dial M for Mobility,
2nd International Workshop on Discrete Algorithms and Methods for13

Mobile Computing and Communications, Dallas, Texas, September 30– October 1 1998.
[12] Y. Lin and S. S. Skiena. Algorithms for square roots of graphs. SIAM
J. Disc. Math., 8:99–118, 1995.
[13] S. T. McCormick. Optimal approximation of sparse Hessians and its
equivalence to a graph coloring problem. Math. Program., 26:153–171,
1983.
[14] V. V. Vazirani. Approximation Algorithms. Springer, 2001. Chapter 5.
Problem
mrng2
mrng2
mrng2
mrng2
mrng2
mrng2

p
1
2
4
6
8
12

χ1
12
12
12
12
12
13

χ3
12
12
12
12
12
13

K
0
1
9
9
9
16

T1
11.5
6.5
4.2
2.8
2.5
1.5

T2
11
6.4
4
2.8
2.5
1.5

T3
0
0
0
0
0
0

Ttot
22.5
12.9
8.3
5.7
5
3

Spar
1
1.7
2.7
4
4.5
7.5

Sseq
0.5
0.9
1.4
2
2.3
3.8

fe144
fe144
fe144
fe144
fe144
fe144

1
2
4
6
8
12

41
40
41
41
43
42

41
40
41
41
43
43

0
3
13
25
20
110

3.8
2.3
1.2
0.8
0.8
0.6

3
2.0
1
0.7
0.7
0.5

0
0
0
0
0
0

6.8
4.3
2.2
1.5
1.5
1.1

1
1.6
3.1
4.5
4.5
6.2

0.6
0.9
1.7
2.5
2.5
3.5

m14b
m14b
m14b
m14b
m14b
m14b

1
2
4
6
8
12

42
43
44
43
44
43

42
43
44
43
46
44

0
0
24
16
28
70

5.6
3.5
1.8
1.7
1.2
0.9

4.8
3.1
1.6
1.4
1
0.8

0
0
0
0
0
0

10.4
6.6
3.4
3.1
2.2
1.8

1
1.6
3.1
3.4
4.7
5.8

0.5
0.9
1.7
1.8
2.6
3.1

ev01
ev01
ev01
ev01
ev01
ev01

1
2
4
6
8
12

3393
3298
3150
3215
2902
2922

3393
3497
3598
3725
3755
3900

0
1514
4555
7786
10113
17872

38.2
18.7
10.3
7.2
5.7
3.9

28.9
14.1
7
5.3
3.9
3

0
7.6
23.4
41.3
52
95

67.1
40.4
40.8
54
61.6
102

1
1.7
1.7
1.2
1.1
0.7

0.6
0.9
0.9
0.7
0.6
0.4

ev02
ev02
ev02
ev02
ev02
ev02

1
2
4
6
8
12

4260
4212
4184
4172
4152
4036

4260
4353
4633
4660
4870
5168

0
1012
4050
6206
8873
16394

144
78.8
46
25
20.6
15.5

111
58.9
28
18.8
15.4
11

0
11.2
45
66
104
188

255
148.9
119
110
140
214.7

1
1.7
2.1
2.3
1.8
1.2

0.6
1
1.2
1.3
1
0.7

Table 2: Performance of Algorithm SD2.

14

Problem
mrng2
mrng2
mrng2
mrng2
mrng2
mrng2

p
1
2
4
6
8
12

χ1
12
12
13
12
12
12

χ2
11
11
10
10
11
10

χ4
11
11
10
10
11
10

K1
0
1
0
6
9
20

K2
0
0
0
0
0
0

T1
9.6
4.8
3.6
2.8
2
1.5

T2
11
5.7
5.4
3.8
2.7
2

T3
8.9
4.7
2.8
2.4
2
1.5

T4
0
0
0
0
0
0

Ttot
29.6
15.3
11.8
9
6.7
5

Spar
1
1.9
2.5
3.3
4.4
6

Sseq
0.7
1.3
1.7
2.3
3.1
4.1

fe144
fe144
fe144
fe144
fe144
fe144

1
2
4
6
8
12

41
41
41
41
41
41

37
38
37
38
37
38

37
38
37
38
37
38

0
0
18
39
55
78

0
0
0
0
0
0

3.2
2.5
1.2
0.9
0.8
0.5

3.8
2.8
1.2
1
1
0.7

2.6
2.1
1
0.7
0.8
0.4

0
0
0
0
0
0

9.6
7.4
3.4
2.6
2.6
1.6

1
1.3
2.8
3.7
3.7
6

0.7
1
2.1
2.7
2.7
4.4

m14b
m14b
m14b
m14b
m14b
m14b

1
2
4
6
8
12

42
43
43
43
44
44

41
41
41
41
41
41

41
41
41
41
41
41

0
0
26
17
31
45

0
0
0
0
0
0

6
3.9
2
1.6
1.2
0.9

8.8
4.9
2.9
1.5
1.3
1

5
3.4
1.7
1
1
0.7

0
0
0
0
0
0

19.8
12.3
6.6
4.1
3.6
2.7

1
1.6
3
4.8
5.5
7.3

0.7
1.2
2.2
3.6
4.1
5.5

ev01
ev01
ev01
ev01
ev01
ev01

1
2
4
6
8
12

3393
3301
3159
3220
2927
2906

3148
3132
3091
3116
3041
3071

3148
3290
3263
3371
3468
3481

0
1687
4171
7913
10932
18264

0
316
863
1475
2675
3040

38.5
18.4
10
7.8
5.3
3.8

44.5
24.7
19
17.5
15.3
17

28.5
13.4
7.3
5.4
4
2.8

0
1.8
4.6
8.5
14.3
16.4

111.5
58.4
41
39.2
39
40.6

1
1.9
2.7
2.8
2.9
2.7

0.7
1.4
2
2.1
2.1
2

ev02
ev02
ev02
ev02
ev02
ev02

1
2
4
6
8
12

4260
4210
4194
4186
4138
4030

4016
4024
4023
4031
4013
4008

4016
4085
4226
4247
4349
4455

0
1019
4154
6196
10732
18392

0
207
644
1354
2007
2819

150
75.2
42.3
25.4
20
14

191
109.2
76.5
48.6
52
55.2

115
57.5
29.2
19
15
10

0
2.7
7.6
16
24
32.7

456
244.7
155.7
109
111
112.5

1
1.9
2.9
4.2
4.1
4.1

0.7
1.4
2.7
3.1
3.1
3

Table 3: Performance of Algorithm ID2.

15

Problem
mrng2
mrng2
mrng2
mrng2
mrng2
mrng2

p
1
2
4
6
8
12

χ1
10
10
10
10
11
11

χ3
10
10
10
10
11
11

K
0
1
6
4
12
10

T1
9.3
5.2
3.3
2.5
2.3
1.7

T2
9
4.8
2.7
2.5
2.3
1.7

T3
0
0
0
0
0
0

Ttot
18.3
10
6
5
4.6
3.4

Spar
1
1.8
3
3.7
4
5.4

Sseq
0.5
0.9
1.6
1.9
2
2.7

fe144
fe144
fe144
fe144
fe144
fe144

1
2
4
6
8
12

35
36
35
36
35
36

35
36
35
36
35
36

0
0
6
15
11
39

4.4
2.8
1.4
1
0.8
0.5

3.7
2.5
1
0.9
0.6
0.5

0
0
0
0
0
0

8.1
5.3
2.4
1.9
1.5
1

1
1.5
3.3
4.3
5.4
8.1

0.5
0.8
1.8
2.3
2.9
4.4

m14b
m14b
m14b
m14b
m14b
m14b

1
2
4
6
8
12

34
37
37
39
37
38

34
37
37
39
37
38

0
0
3
7
10
16

5.4
2.9
1.9
1.6
1.2
1

4.6
2.6
1.6
1.3
1.1
0.8

0
0
0
0
0
0

10
5.5
3.5
2.9
2.3
1.8

1
1.8
2.9
3.5
4.3
5.6

0.5
1
1.5
1.9
2.3
3

ev01
ev01
ev01
ev01
ev01
ev01

1
2
4
6
8
12

2148
1969
1915
2429
1822
2103

2148
2012
2070
2618
2286
2687

0
135
366
538
1031
1424

35.5
19.2
9.5
7.7
5.6
4

29
15.4
7.6
6
4.3
3

0
0.7
1.9
2.8
6.1
7.5

64.5
35.3
19
16.5
16
14.5

1
1.8
3.4
3.9
4
4.4

0.6
1
1.9
2.2
2.2
2.5

ev02
ev02
ev02
ev02
ev02
ev02

1
2
4
6
8
12

2697
2626
2590
2584
2626
2536

2697
2653
2736
2818
2967
3049

0
83
340
595
960
1569

134.4
74.4
40
29
20
15.7

112.6
60.1
32.2
21
15.6
11.4

0
1
4.4
6.7
10
18.6

247
135.5
76.5
56.7
45.5
45.7

1
1.8
3.2
4.4
5.4
5.4

0.5
1
1.8
2.4
3
3

Table 4: Performance of Algorithm SD 32 .

Problem
ev01
ev01
ev01
ev01
ev01
ev01

p
1
2
4
6
8
12

χ1
4077
4084
4043
4132
3998
4019

χ3
4077
4088
4051
4144
4015
4042

K
0
59
139
572
388
587

T1
38
20
10
7
6
4

T2
28
15
7.8
5.5
4.3
2.8

T3
0
0.3
0.8
3.5
2
3.2

Ttot
66
35.2
19
16
12
10

Spar
1
1.9
3.5
4.1
5.5
6.6

Sseq
0.6
1.1
2
2.4
3.2
3.8

ev02
ev02
ev02
ev02
ev02
ev02

1
2
4
6
8
12

5564
5581
5518
5553
5576
5514

5564
5584
5525
5562
5583
5536

0
43
187
273
558
522

148.3
82
41
26
21
15

109
60
31
19
16
10

0
0.6
2.3
3.3
7.4
6

257.3
143
75
49
44
31

1
1.8
3.4
5.2
5.8
8.3

0.6
1
2
3
3.4
4.8

Table 5: Performance of Algorithm RSD2, q = 1/20.

16

Problem
ev01
ev01
ev01
ev01
ev01
ev01

p
1
2
4
6
8
12

χ1
3719
3706
3648
3736
3601
3628

χ2
3169
3161
3183
3205
3167
3195

χ4
3169
3173
3217
3326
3236
3291

K1
0
171
483
998
1272
1857

K2
0
21
152
193
342
447

T1
36.6
19
10
8
5.7
3.7

T2
42.5
27
21
21
19
18.7

T3
27
15
7.5
5.3
4
2.6

T4
0
0.1
0.8
1
2
2.2

Ttot
106.1
61
40
35
30
27

Spar
1
1.7
2.7
3
3.5
3.9

Sseq
0.7
1.3
2
2.3
2.6
2.9

ev02
ev02
ev02
ev02
ev02
ev02

1
2
4
6
8
12

4919
4871
4860
4865
4869
4839

4024
4027
4025
4040
4032
4031

4024
4066
4082
4084
4100
4104

0
208
526
787
1577
1953

0
43
127
138
207
320

132.7
72
37
27
22
15

150
98
62
59
61
75

97.3
52
27
19
15
12

0
0.5
1.5
1.7
2.3
4

380
223
128
108
101
106

1
1.7
3
3.5
3.8
3.6

0.7
1.3
2.2
2.6
2.8
2.7

Table 6: Performance of Algorithm RID2, q = 1/6.

Problem
ev01
ev01
ev01
ev01
ev01
ev01

p
1
2
4
6
8
12

χ1
2392
2242
2189
2582
2170
2496

χ3
2392
2242
2190
2583
2225
2528

K
0
23
46
55
138
191

T1
37
18
9.7
7.6
5.8
3.8

T2
30
15
7.8
5.4
4.8
3

T3
0
0.1
0.3
0.3
0.7
1.3

Ttot
67
33
18
13
11
8

Spar
1
2
3.7
5.1
6.1
8.4

Sseq
0.6
1.1
2
2.8
3.4
4.6

ev02
ev02
ev02
ev02
ev02
ev02

1
2
4
6
8
12

3142
3049
3076
3018
3106
2995

3142
3049
3080
3060
3125
3081

0
18
51
123
155
214

140
71
38
25
20
15

117
56
31
21
16
12

0
0.2
0.6
1.4
1.7
2.4

257
127
69
48
37
29

1
2
3.7
5.3
6.9
8.8

0.5
1.1
2
2.9
3.8
4.8

Table 7: Performance of Algorithm RSD 32 , q = 1/6.

17

Graph Coloring on Coarse Grained Multicomputers∗
Assefaw Hadish Gebremedhin†Jens Gustedt§

Isabelle Guérin Lassous‡
Jan Arne Telle

Abstract
We present an efficient and scalable Coarse Grained Multicomputer(CGM) coloring algorithm that colors a graph G with at most ∆ + 1
colors where ∆ is the maximum degree in G. This algorithm is given
in two variants: randomized and deterministic. We show that on ap-processor CGM model the proposed algorithms require a paralleltime of O( |E|p ) and a total work and overall communication cost ofO(|E|).These bounds correspond to the average case for the randomized version and to the worst case for the deterministic variant.
key words: graph algorithms, parallel algorithms, graph coloring,
Coarse Grained Multicomputers.

∗

Research supported in part by The Aurora Programme, a France-Norway Collaboration Research Project of The Research Council of Norway, The French Ministry of Foreign
Affairs and The Ministry of Education, Research and Technology.
†
Department of Informatics, University of Bergen, 5020 Bergen, Norway. {assefaw,
telle}@ii.uib.no‡
INRIA Rocquencourt, France. Isabelle.Guerin-Lassous@inria.fr§
LORIA & INRIA Lorraine, France. Jens.Gustedt@loria.fr

1

1

Introduction

The graph coloring problem deals with the assignment of positive integers(colors) to the vertices of a graph such that adjacent vertices do not get thesame color and the number of colors used is minimized. A wide range of realworld problems, among others, timetabling and scheduling [20], frequency
assignment [8], register allocation [2], and efficient estimation of sparse matrices in optimization [4], have successfully been modeled using the graphcoloring problem. Besides modeling real world problems, graph coloring
plays a crucial role in the field of parallel computation. In particular, whena computational task is modeled using a graph where the vertices represent the subtasks and the edges correspond to the relationship among them,graph coloring is used in dividing the subtasks into independent sets thatcan be performed concurrently.
The graph coloring problem is known to be NP-hard [9], making heuristic approaches inevitable in practice. There exist a number of sequential
graph coloring heuristics that are quite effective in coloring graphs encountered in practical applications. See [4] for some of the popular heuristics.
However, due to their inherent sequential nature, these heuristics are difficult to parallelize. In fact, coloring the vertices of a graph in a given
order where each vertex is assigned the smallest color that has not been
given to any of its neighbors is shown to be P-complete [12]. Consequently,
parallel graph coloring heuristics different from the effective sequential coloring heuristics had to be suggested. One of the important contributionsin this regard is the parallel maximal independent set finding algorithmof Luby [21] and the coloring algorithm based on it. Subsequently, Jonesand Plassmann [16] improved Luby’s algorithm and in addition used graphpartitioning as a means to achieve a distributed memory coloring heuristic
based on explicit message-passing. Unfortunately, Jones and Plassmann did
not get any speedup from their experimental studies. Later, Allwrignt etal. [1] performed a comparative study of the implementations of the JonesPlassmann algorithm and a few other variations and reported that none of
the algorithms included in the study yielded any speedup. The justification
for the use of these parallel coloring heuristics has been the fact that they
enabled solving large-scale problems that could not otherwise fit onto thememory of a sequential machine.
Despite these discouraging experiences, Gebremedhin and Manne [10]recently proposed a shared memory parallel coloring algorithm that yields
good speedup. Their theoretical analysis using the PRAM model showsthat the algorithm is expected to provide an almost linear speedup and2

experimental results conducted on the Origin 2000 supercomputer usinggraphs that arise from finite element methods and eigenvalue computations
validate the theoretical analysis.
The purpose of this paper is to make this successful approach feasible
for a larger variety of architectures by extending it to the Coarse Grained
Multicomputer (CGM) model of parallel computation [6]. The CGM model
makes an abstraction of the interconnection network among the processorsof a parallel computer (or network of computers) and captures the efficiency
of a parallel algorithm using only a few parameters. Several experiments
show that the CGM model is of practical relevance: implementations ofalgorithms formulated in the CGM model in general turn out to be portable,
predictable, and efficient [13, 14].
In this paper we propose a CGM coloring algorithm that colors a graph
G with at most ∆ + 1 colors where ∆ is the maximum degree in G. The
algorithm is given in two variants: one randomized and the other deterministic. We show that the proposed algorithms require a parallel time ofO( |E|p ) and a total work and overall communication cost of O(|E|). Thesebounds correspond to the average case for the randomized version and to
the worst case for the deterministic variant.
The remainder of this paper is organized as follows. In Section 2 wereview the CGM model of parallel computation and the graph coloring
problem. In Section 3 we discuss a good data organization for our CGM
algorithms and present the randomized variant of the algorithm along with
its various subroutines. In Section 4 we provide an average case analysisof the randomized algorithm’s time and work complexity. In Section 5 weshow how to de-randomize our algorithm to achieve the same good time andwork complexity also in the worst case. Finally, in Section 6 we give someconcluding remarks.

2

Background2.1

Coarse grained models of parallel computation

In the last decade several efforts have been made to define models of parallel (or distributed) computation that are more realistic than the classical
PRAM models; see [7] or [19] for an overview of PRAM models. In contrast
to the PRAM models that suppose that the number of processors p is polynomial in the input size N , the new models are coarse grained, i.e., they
assume that p and N are orders of magnitude apart. Due to this assumption, the coarse grained models map much better on existing architectures3

where in general the number of processors is in the order of hundreds andthe size of the data to be handled could be in the order of billions.
The introduction of the Bulk Synchronous Parallel (BSP) bridging model
for parallel computation by Valiant [24] marked the beginning of the increasing research interest in coarse grained parallel computation. The BSP model
was later modified along different directions. For example, Culler et al. [5]suggested the LogP model as an extension of Valiant’s BSP model in whichasynchronous execution was modeled and a parameter was added to betteraccount for communication overhead. In an effort to define a parallel computation model that retains the advantages of coarse grained models whileat the same time is simple to use (involves few parameters), Dehne et al. [6]suggested the CGM model.
The CGM model considered in this paper is well suited for the design ofalgorithms that are not too dependent on a particular architecture and ourbasic assumptions of the model are listed below.
• The model consists of p processors and all the processors have thesame size M = O(N/p) of memory, where N is the input size.
• An algorithm on this model proceeds in so-called supersteps. A superstep consists of one phase of local computation and one phase ofinterprocessor communication.
• The communication network between the processors can be arbitrary.
The goal when designing an algorithm in this model is to keep the sum
total of the computational cost per processor, the overall communication
cost, and idle time of each processor within T /s(p), where T is the runtimeof the best sequential algorithm on the same input, and the speedup s(p) isa function that should be as close to p as possible.
To achieve this, it is desirable to keep the number of supersteps of suchan algorithm as low as possible, preferably within O(M ). The rationale
here lies in the fact that, among others, the message startup-cost and the
bandwidth of an architecture determine the communication overhead. In
each superstep, a processor may need to do at most O(p) communications
and hence a number of supersteps of O(M ) ensures that the total startupcost is at most O(M p) = O(N ) and therefore lies within the complexity
bound of the overall computational cost we anticipate for such an algorithm.
The bandwidth restriction of a specific platform must still be observed,
and here the best strategy is to reduce the communication volume as muchas possible. See [13] for an overview of algorithms, implementations andexperiments on the CGM model.
As a legacy from the PRAM model, it is usually assumed that the num4

ber of supersteps should be polylogarithmic in p. However, the assumption
seems to have no practical justification. In fact, there is no known relationship between the coarse grained models and the complexity classes N C k .
In practice, algorithms that simply ensure a number of supersteps that is afunction of p (but not of N ) perform quite well [11].
To be able to organize the supersteps well, we assume that each processor
can store a vector of size p for every other processor. Thus, the following
inequality is assumed throughout this paper,p2 < M.

2.2(1)

Graph coloring

A graph coloring is a labeling of the vertices of a graph G = (V, E) withpositive integers, called colors, such that adjacent vertices do not obtain thesame color. It can equivalently be viewed as searching for a partition of
the vertex set of the graph into independent sets. The primary objectivein the graph coloring problem is to minimize the number of colors used.
Even though coloring a graph with the fewest number of colors is an NPhard problem, in many applications coloring using a bounded number ofcolors, possibly far from the minimum, may suffice. Particularly in many
parallel graph algorithms, a bounded coloring (partition into independentsets) is needed as a subroutine. For example, graph coloring is used in thedevelopment of a parallel algorithm for computing the eigenvalues of certain
matrices [22] and in parallel partial differential equation solvers [1].
One of the simplest and yet quite effective sequential heuristics for graphcoloring is the greedy algorithm that visits the vertices of the graph in someorder and in each visit assigns a vertex the smallest color that has not beenused by any of the vertex’s neighbors. It is easy to see that, for a graph
G = (V, E), such a greedy algorithm always uses at most ∆+1 colors, where∆ = maxv∈V {degree of v}. In Greenlaw et al. [12], a restricted variant of
the greedy algorithm in which the ordering of the vertices is predefined, andthe algorithm is required to respect the given order, is termed as Lexicographically First ∆ + 1-coloring (LF∆ + 1-coloring). We refer to the casewhere this restriction is absent and where the only requirement is that theresulting coloring uses at most ∆ + 1 colors, simply as ∆ + 1-coloring.
LF∆ + 1-coloring is known to be P-complete [12]. But for special classesof graphs, some N C algorithms have been developed for it. For example,
Chelbus et al. [3] show that for tree structured graphs LF∆ + 1-coloringis in N C. In the absence of the lexicographically first requirement, a few
N C algorithms for general graphs have been proposed. Luby [21] has given5

an N C algorithm for ∆ + 1-coloring by reducing the coloring problem tothe maximal independent set problem. Moreover, Karchmer and Naor [17],Karloff [18], and Hajnal and Szemerédi [15] have each presented different
N C algorithms for Brook’s coloring (a coloring that uses at most ∆ colors
for a graph whose chromatic number is bounded by ∆). Earlier, Naor [23]had established that coloring planar graphs using five colors is in N C.
However, all of these N C coloring algorithms are mainly of theoreticalinterest as they require a polynomial number of processors, whereas, in reality, one has only a limited number of processors on a given parallel computer.
In this regard, Gebremedhin and Manne [10] have recently shown a practical and effective shared memory parallel ∆ + 1-coloring algorithm. Theyshow that distributing the vertices of a graph evenly among the available
processors and coloring the vertices on each processor concurrently, while
checking for color compatibility with already colored neighbors, creates very
few conflicts. More specifically, the probability that a pair of adjacent vertices are colored at exactly the same instance of the computation is quite
small. On a somewhat simplified level, the algorithm of Gebremedhin andManne works by tackling the list of vertices numbered from 1 to n in a‘round robin’ manner. At a given time t, where 1 ≤ t ≤ r and r = d np e,
processor Pi colors vertex (i − 1) · r + t. The shared memory assumptions
ensure that Pi may access the color information of any vertex at unit costof time. Adjacent vertices that are in fact handled at exactly the same time
are the only causes for concern as they may result in conflicts. Gebremedhinand Manne show that the number of such conflicts is small on expectation,and that conflicts can easily be resolved a posteriori. Their analysis of theresulting algorithm using the PRAM model shows that the algorithm colors
a general graph G = (V, E) with ∆ + 1 colors in expected time O(|E|/p)pwhen the number of processors p is such that p ≤ |V |/ 2|E|.
However, in a distributed memory setting, the most common case in ourtarget model CGM, one has to be more careful about access to data locatedon other processors.

3

A CGM ∆ + 1-coloring algorithm

We start this section by discussing how we distribute the input graph among
the available processors for our CGM ∆ + 1-coloring algorithms. Then, therandomized variant of our algorithm is presented in a top-down fashion,
starting with an overview and filling the details as the presentation proceeds.

6

3.1

Data distribution

In general, a good data organization is crucial for the efficiency of a distributed memory parallel algorithm. For our CGM-coloring algorithm inparticular, the input graph G = (V, E) is organized in the following manner.
• Each processor Pi (1 ≤ i ≤ p) is responsible for a subset U i of theSvertices (V = pi=1 Ui ). With a slight abuse of notation, the processor
hosting a vertex v is denoted by Pv .
• Each edge e = {v, w} ∈ E is represented as arcs (v, w) stored at P v ,and (w, v) stored at Pw .
• For each arc (v, w) processor Pv stores the identity of Pw and thus thelocation of the arc (w, v). This is to avoid a logarithmic blow-up dueto the search for Pw .
• The arcs are sorted lexicographically and stored as a linked list pervertex.
In this data distribution, we require that the degree of each vertex be lessthan D = d Np e, where N = |E|. Vertices with degree greater than D aretreated in a separate preprocessing step.
If the input of the algorithm is not of the desired form, it can be efficiently
transformed into one by carrying out the following steps.
• Generate two arcs for each edge as described above,• Radix sort (see [13] for a CGM radix sort) the list of arcs such thateach processor receives the arc (v, w) if it is responsible for vertex w,• Let every processor note its identity on these arcs,• Radix sort the list of arcs such that every processor receives its properarcs (arc (v, w) if it is responsible for vertex v).

3.2

The algorithm

As the time complexity of sequential ∆ + 1-coloring is linear in the numberof edges |E|, our aim is to design a parallel algorithm in CGM with O( |E|p )work per processor and O(|E|) overall communication cost. In an overview,
our CGM coloring algorithm consists of two phases, an initial and a mainrecursive phase; see Algorithm 1.
In the initial phase, the subgraph induced by the vertices with degree
greater than d Np e is colored sequentially on one of the processors. Clearly,
there are at most p such vertices since otherwise we would have more thanN edges in total. Thus the subgraph induced by these vertices has at mostp2 edges. Since p2 is assumed to be less than M , the induced subgraph fits
7

Algorithm 1: ∆ + 1-coloring on a CGM with p processors
Input: Base graph G = (V, E), the subgraph H induced by verticesof degree greater than D = dN/pe, the lists F v of forbidden
colors of vertices v ∈ V .
Output: A valid coloring of G = (V, E) with at most ∆ + 1 colors.
initial phase Sequential∆ + 1Coloring(H,{F v }v ) (see Algorithm 3);
main phase ParallelRecursive∆ + 1Coloring(G,{F v }v ) (see Algorithm 2);
P1
P2
P3
P4
P5
P6
t=1

t=2

t=3

t=4

Figure 1: Graph on 72 vertices distributed onto 6 processors and 4 timeslots.
on a single processor (say P1 ) and a call to Algorithm 3 colors it sequentially.
Algorithm 3 is also used in a situation other than coloring such vertices. We
defer the discussion on the details of Algorithm 3 to Section 3.2.1 where thesituation that calls for its second use is presented.
The main part of Algorithm 1 is the call to Algorithm 2 which recursively
colors any graph G such that the maximum degree ∆ ≤ M . The basic ideaof the algorithm is based on placing the vertices residing on each processor
into different timeslots. The assignment of timeslots to the vertices gives
rise to two categories of edges. The first category consists of edges whichconnect vertices having the same timeslot. We call these edges bad and allother edges good. Figure 1 shows an example of a graph distributed on 6
processors and 4 timeslots in which bad edges are solid and good edges aredashed.
In a nutshell, Algorithm 2 proceeds timeslot by timeslot where in eachtimeslot the graph defined by the bad edges and the vertices incident onthem is identified and the algorithm is called recursively with the identified
graph as input while the rest of the input graph is colored concurrently.
In Algorithm 2, while partitioning the vertices into k timeslots, where8

Algorithm 2: Parallel Recursive ∆ + 1-coloring
Input: Subgraph G0 = (V 0 , E 0 ) of a base graph G = (V, E) with
M 0 edges per processor such that ∆G0 ≤ M 0 , M the initial
input size per processor, lists Fv of forbidden colors for thevertices.
Output: A valid coloring of G0 with at most ∆G + 1 colors.
Nbase case if (|G0 | < kp2 ) thenSequential∆ + 1Coloring(G0 ,{Fv }v ) (see Algo-rithm 3);else
high degree
HandleHighDegreeVertices(G0 ,{Fv }v ,2k) (see Algorithm 5);group vertices
foreach Pi do
Let Ui,t for t = 1, . . . , k be result of the callGroupVerticesIntoTimeslots(V 0 ,k), (see Algorithm 6);
For each vertex v denote the index of its timeslot by t v ;foreach arc (v, w) do collect the timeslot t v in a send buffer
for Pw ;Send out the tuples (w,tv );Receive the timeslots from the other processors;
for t = 1 to k doforeach processor Pi do
identify conflicts
Consider all arcs e = (v, w) with v ∈ U i,t and tv = tw = t;Name this set Si and consider the vertices VSi that havesuch an arc;balance load
recursecolor vertex
send messages

receive messages

Grec =Balance( pi=1 VSi , pi=1 Si ) (see Algorithm 7);
ParallelRecursive∆ + 1Coloring(G rec ,{Fv }v );foreach processor Pi do
foreach uncolored vertex v with t v = t do Color v withleast legal color;foreach arc (v, w) with v ∈ Ui , tv = t and tw > t do
Collect the color of v in a send buffer for P w ;Send out the tuples (w,color of v);
Receive the colors from the other processors;foreach received tuples (w, color of v) do add color of v
to Fw ;SS

9

1 < k ≤ p, we would like to achieve as even a distribution as possible. The
call to Algorithm 6 in line group vertices does this by using the degreeof each vertex as a criterion. This randomized algorithm is presented in
Section 3.2.2 where the issue of load balancing is briefly discussed. Prior tocalling Algorithm 6, vertices with ‘high degrees’ that would otherwise resultin an uneven load balance are treated separately; see line high degree.
The algorithm for treating high degree vertices, Algorithm 5, is presented
in Section 3.2.2.
Notice that an attempt to concurrently color vertices incident on a badedge may result in an inconsistent coloring (conflict). In a similar situation,
Gebremedhin and Manne, in their shared memory formulation, tentatively
allow such conflicts and resolve eventual conflicts in a later sequential phase.
The success of their approach lies in the fact that the expected size of theedges in conflict is relatively small. In our case, we deal with the potential
conflicts a priori. We first identify the subgraphs that could result in conflict and then color these subgraphs in parallel recursively until their unionis small enough to fit onto the memory of a single processor. See lines identify conflicts and recurse in Algorithm 2. Note that, in general, someprocessors may receive more vertices than others. We must ensure thatthese recursive calls do not produce a blow-up in computation and communication. In order to ensure that the subgraph that goes into recursion isevenly distributed among the processors, a call to Algorithm 7 is made atline balance load. Algorithm 7 is discussed in Section 3.2.2.
In the recursive calls one must handle the restrictions that are imposedby previously colored vertices. We extend the problem specification andassume that a vertex v also has a list F v of forbidden colors that initiallyis empty. An important issue for the complexity bounds is that a forbidden
color is added to Fv only when the knowledge about it arrives on P v . The
list Fv as a whole will only be touched once, namely when v is finally colored.
Observe also that the recursive calls in line recurse need not be synchronized. In other words, it is not necessary (nor desired) that the processors
start recursion at exactly the same moment in time. During recursion, when
the calls reach the communication phase of the algorithm, they will automatically be synchronized in waiting for data from each other.
Clearly, the subgraph defined by the good edges and their incident vertices can be colored concurrently by the available processors. In particular,
each processor is responsible for coloring its set of vertices as shown in line
color vertex of Algorithm 2. In determining the least available color to avertex, each processor maintains a Boolean vector Bcolors. This vector isindexed with the colors and initialized with all values set to “true”. Then
10

Algorithm 3: Sequential ∆ + 1Coloring
Input: M the initial input size per processor, subgraph G 0 = (V 0 , E 0 )of a base graph G = (V, E) with |E 0 | ≤ M and lists Fv offorbidden colors for the vertices.
find allowed foreach processor Pi do
Let Ui0 = Ui ∩ V 0 be the vertices of G0 that are stored on Pi ;For each v ∈ Ui0 let d(v) be the degree of v in G0 ;Av = ComputeAllowedColors(v,d(v),{F v }v ) (see Algorithm 4);
Communicate E 0 and all lists Av to P1 ;color sequentially for processor P1 do
Collect the graph G0 together with the lists Av ;Color each vertex in G0 with least available color;
Send the resulting colors back to the corresponding processors;communicate foreach processor Pi do
Inform all neighbors of Ui0 of the colors that have been assigned;
Receive the colors from the other processors and update the lists
Fv accordingly;

when processing a vertex v, the entries of Bcolors corresponding to v’s listof forbidden colors are set to “false”. After that, the first item in Bcolors
that still is true is looked for and chosen as the color of v. Then, the vector
is reset by assigning all its modified values the value “true” for future use.
After a processor has colored a vertex, it communicates the color information to processors hosting a neighbor. In each timeslot the messages tothe other processors are grouped together, see send messages and receive
messages. This way at most p − 1 messages are sent per processor pertimeslot.
3.2.1

The base case

The base case of the recursion is handled by a call to Algorithm 3 (see basecase in Algorithm 2). Note that the sizes of the lists F v of forbidden colors
that the vertices might have collected during higher levels of recursion may
actually be too large and their union might not fit on a single processor.
To handle this situation properly, we proceed in three steps as shown in
Algorithm 3. Notice that Algorithm 3 is the same routine called in theinitial phase of Algorithm 1.
In the step find allowed, for each vertex v ∈ V 0 a short list of allowed
colors Av is computed. Observe that a vertex v can always be colored using
one color from the set {1, 2, . . . , d(v)+1}, where d(v) is the degree of v. Hencea list of d(v) + 1 allowed colors suffices to take all restrictions of forbidden
11

Algorithm 4: Compute Allowed Colors
Input: v together with its actual degree d(v) and its (unordered)
list Fv of forbidden colors; A Boolean vector colors with allvalues set to true.
Output: a sorted list Av of the least d(v) + 1 allowed colors for v
foreach c ∈ Fv do Set colors[c] = f alse;for (c = 1; |Av | < d(v); + + c) do if colors[c] then Av = Av + c;
foreach c ∈ Fv do Set colors[c] = true;
colors into account. Using a similar technique as described in color vertexof Algorithm 2, we can obtain a sorted list A v of allowed colors for v in time
proportional to |Fv | + d(v). This is done by the call to Algorithm 4 in line
find allowed. Then in the step color sequentially, the vertices of the
input graph are colored sequentially using their computed lists of allowed
colors. In the final step communicate, the color information of the vertices
is communicated.
3.2.2

Load balancing

In this section we address the issue of load balancing. In Algorithm 2,three matters that potentially result in an uneven load balance are (i) highvariation in the degrees of the vertices, (ii) high variation in the sum of thedegrees in the timeslots, and (iii) the recursive calls on the subgraphs that
go into recursion. The following three paragraphs are concerned with thesepoints.
Handling high degree vertices Whereas for the shared memory algorithm differences in degrees of the vertices that are colored in parallel justcauses a slight asynchrony in the execution of the algorithm, in a CGM setting it might result in a severe load imbalance and even in memory overflow
of a processor.
Line group vertices of Algorithm 2 groups the vertices into k ≤ ptimeslots of about equal degree sum. If the variation in the degrees of thevertices is too large, such a grouping would not be even. For example, if wehave one vertex of very large degree, it would always dominate the degree
sum of its time slot thereby creating imbalance. So, we have to make surethat the degree of each vertex is fairly small, namely smaller than dM 0 /qe
where q is a parameter of Algorithm 5. Observe that the notion of ‘small’degree depends on the input size M 0 and thus may change during the courseof the algorithm. This is why we need to have the line high degree in every

12

Algorithm 5: Handle High Degree Vertices
Input: Subgraph G0 = (V 0 , E 0 ) of a base graph G = (V, E) with M 0edges per processor such that ∆G0 ≤ M 0 , lists Fv of forbidden
colors for the vertices and a parameter q.
foreach processor Pi do
find all v ∈ Ui with degree higher than M 0 /q (Note: all degrees
are less than N/p);send the names and the degrees of these vertices to P 1 ;for processor P1 do
Receive lists of high degree vertices;Group these vertices into k 0 ≤ q timeslots W1 , ..., Wk0 of at mostp vertices each and of a degree sum of at most 2N/p for eachtimeslot;
Communicate the timeslots to the other processors;foreach processor Pi do
Receive the timeslots for the high degree vertices in U i ;Communicate these values to all the neighbors of these vertices;
Receive the corresponding information from the other processors;
Compute Et,i for t = 1, . . . , k 0 where one endpoint is in Wt ∪ Ui ;for t = 1 to k 0 do
S
Let Et = 1≤i≤p Et,i and denote by Gt = (Wt , Et ) the induced
subgraph of high degree vertices of timeslot t;
Sequential∆ + 1Coloring(Gt ,{Fv }v ) (see Algorithm 3);

recursive call and not only at the top level call. Note that q is a multiple ofk, the number of timeslots of Algorithm 2.
Thus, the high degree vertices that we indeed have to treat in eachrecursive call are those vertices v with dM 0 /qe < deg(v) ≤ M 0 . Such vertices
are handled using Algorithm 5, which essentially divides the set of highdegree vertices into k 0 ≤ q timeslots and colors each of the subgraphs inducedby these timeslots sequentially.
Grouping vertices into timeslots Algorithm 6 partitions the vertices
into k timeslots. It does so by first dividing the set of vertices into groupsof size k and then distributing the vertices of each group into the distinct
timeslots. Observe that no communication is required during the course ofthis algorithm.
The partition obtained with this algorithm is relatively balanced.
Lemma 1 On each processor P , the difference of the degree sums of thevertices in any two timeslots is at most the maximum degree over all vertices
13

Algorithm 6: Group Vertices Randomly into Timeslots
Input: V 0 the set of vertices, kforeach processor Pi do
Radix sort its vertices according to their descending degrees;
Let v1 , . . . , vs be this order of the vertices;
for i = 0, . . . , d ks e − 1 do
Let j1 , . . . , jk be a random permutation of the values 1, . . . , k;
Assign vik+1 , . . . , v(i+1)k to timeslots j1 , . . . , jk respectively;that P holds.
Proof: Since the vertices are considered in descending order of their degrees, the difference in degree sums between two timeslots is maximized
when one of the timeslots always receives the vertex with the highest degreein the group and the other the smallest. In group i, the vertex of highest degree is vik+1 and the one of smallest degree is v(i+1)k . Thus we can estimate
the difference as follows:d ks e−1

d ks e−1

d ks e−1

X

X

X

i=0

deg(vik+1 )−deg(v(i+1)k ) ≤

i=0

i=0

d ks e−2deg(vik+1 )−Xdeg(v(i+1)k+1 )i=0

which is in turn bounded by deg(v1 ), where v1 has the maximum degree
over all vertices that P holds.
From Lemma 1 and from the fact that we do not have high degree vertices, it follows that the sum of the degrees of the vertices in any timeslot is0
3M 0
between M
2k and 2k .
Balancing during recursion In Algorithm 2, unless proper attention ispaid, the edges of the subgraph that goes into recursion may not be evenly
distributed among the processors. To address this, we suggest an algorithm
that ensures that Grec , the graph that goes into recursion in Algorithm 2,
is evenly distributed among the processors. See Algorithm 7.
Obviously Algorithm 7 runs in time proportional to the input size oneach processor and has a constant number of supersteps.

4

Average case analysis

In this section we provide an average case analysis of Algorithm 2. In
Section 5 we show how to replace the randomized algorithm, Algorithm 6,

14

Algorithm 7: Balance
Input: Graph G0 = (V 0 , E 0 ), such that each v ∈ V 0 has degG0 (v) ≤|E 0 |/p.
Output: A redistribution of V 0 and E 0 on the processors such thateach processor handles no more than M 0 = 2|E 0 |/p edges.
Initialize a distributed array Deg indexed by V 0 that holds the degrees of all vertices;
Do a prefix sum on Deg and store this sum in a similar array P re;
foreach processor Pi do
foreach v ∈ V 0 ∩ Ui do
Let j ∈ {1, ..., p} be such that jM 0 ≤ P re[v] < (j + 1)M 0 ;Send v and its adjacent edges to processor P j ;foreach processor Pi do
receive the corresponding vertices and edges

by a deterministic one. Every lemma in this section refers to Algorithm 2unless stated otherwise.
Lemma 2 For any edge {v, w}, the probability that t v = tw is at most k1 .
Proof: Consider Algorithm 6. We distinguish between two cases. The firstis the case where v and w reside on different processors. In this case, thechoices for the timeslots of v and w are clearly independent, implying thatthe probability that w is in the same timeslot as v is k1 .
The same argument applies for the case where v and w reside on the sameprocessor but are not processed in the same group. Whenever they are in
the same group, they are never placed into the same timeslot. Therefore,
the overall probability is bounded by k1 .
Lemma 3 The expected sum total of the number of edges of all subgraphs0
going into recursion in recurse is at most |Ek | .
Proof: The expected total number of edges going into recursion is equal
to the expected total number of bad edges. The latter is in turn equal toP|E 0 |e∈E 0 prob(e is bad), which by Lemma 2 can be bounded by k .
Lemma 4 The expected overall size of the subgraphs at the ith recursion
level is at most N/k i , with at most M/k i per processor.
Proof: Notice that the choices of timeslots between two successive recursion
levels may not be independent. However, the dependency that may occur
actually reduces the number of bad edges even more. This can be seen froma similar argument as that of Lemma 2: vertices that are in the same groupof the degree sequence in Algorithm 6 are forced to be separated into two15different timeslots. For all others, the choices are again independent.
Thus, the total number of edges going into recursion can be immediately
bounded by N/k i . The fact that it is also balanced across the processors isdue to Algorithm 7.
Lemma 5 The expected sum total of the sizes of all the subgraphs handled
by any processor during Algorithm 2 (including all recursions) is O(M ).
Proof: By Lemma 4, the expected sum of the sizes of these graphs isbounded by∞Xkk −i M =M ≤ 2M,(2)
k−1
i=0
for all k ≥ 2. Thus, the total expected size in all the steps per processor isO(M ).
Lemma 6 For any 1 < k ≤ p, the expected number of supersteps is at mostquadratic in p.
Proof: The expected recursion depth of our algorithm is the minimum
value d such that N/k d ≤ M = N/p, which implies k d ≥ p, i.e. d = dlog k pe.
The total number of supersteps in each call (including the supersteps in
Algorithm 5) is c · k, for some constant c ≥ 1. The constant c captures thefollowing supersteps:• some to handle high degree vertices,• one to propagate the chosen timeslots,• some to balance the edges inside each timeslot,• one to propagate the colors for each timeslot.
Thus, the total number of supersteps on recursion level i is c · k i and the
expected number of supersteps is bounded as follows
dlog k peXc · k i ≤ c · k logk p+1 = c · k · p.

(3)

i=1

Lemma 7 The expected overall work involved in base case is O(M ).
Proof: Algorithm 3 on input G0 = (V 0 , E 0 ) and lists of forbidden colors Fv
has overall work and communication cost proportional to |G 0 | and the sizeof the lists Fv .
There are k dlogk pe expected calls to Algorithm 3 in Algorithm 2; thereforeNdlog k pe N ≤ k 1+logk p N ≤P1 is expected to handle k dlogk pe kp
2 edges and k
kp2
kp2
N
kp kp2 = M . This implies an expected work and communication cost ofO(M ) for base case.
16

Lemma 8 The expected overall work per processor involved in high degree
is O(M ).
Proof: In Algorithm 2, in the first call to Algorithm 5 (M 0 = M ), every
processor holds at most q = 2k high degree vertices (i.e vertices v of degreeNdegv (G) such that pq< degv (G) ≤ Np ). Otherwise, it would hold more than(M/q) · q = M edges. So, overall, there are at most p · q such vertices for the
first level of recursion. Processor P 1 distributes these O(p2 ) vertices onto k 0
timeslots such that each timeslot has a degree sum of at most 2N/p = 2M .
Thus, each timeslot induces a graph of expected size 2M/k 0 . Subsequently,
sequential ∆+1-coloring is called for the subgraph induced by each timeslot,
for total work O(M ) = O(M 0 ).
By induction we see that in the ith level of recursion, if a vertex v is of0high degree, its degree degv (Grec ) has to be Mq < degv (Grec ) ≤ M 0 . Using
the same argument as the one above, it can be shown that the total work
to handle these vertices is O(M 0 ).
From Lemma 5, the total expected work in all the steps per processor isO(M ).
Lemma 9 The expected overall work per processor involved ingroup vertices is O(M ).
Proof: Observe that the radix sort can be done in O(M 0 ), since the sortkeys are less than M 0 . The random permutations can easily be computed
locally in linear time.
Again, Lemma 5 proves the claim.
Theorem 1 For any 1 < k ≤ p, the expected work, communication, andidle time per processor of Algorithm 2 is within O(M ). In particular, theexpected total run-time per processor is O(M ).
Proof: From Lemma 6 we see that the expected number of supersteps isO(p2 ); hence by inequality (1) the expected communication overhead generated in all the supersteps is O(M ).
We proceed by showing that the work and communication that a processor has to perform in Algorithm 2 is a function of the number of edges
on that processor, i.e. M . Inserting a new forbidden color into an unsorted
list Fv can be done in constant time. Since an edge contributes an item tothe list of forbidden colors of one of its incident vertices at most once, thesize of such a list is bounded by the degree of the vertex. Thus, the total
size of these lists on any of the processors will never exceed the input size
M 0 (recall that vertices of degree greater than Np have been handled in the
preprocessing step).

17

As discussed in Section 3.2, a Boolean vector Bcolors is used in determining the color to be assigned to a vertex. In the absence of high degree0vertices no list Fv will be longer than Mq and hence the size of Bcolors
0need not exceed Mq + 1. Even when this restriction is relaxed, as shown in
Section 3.2.2, we need at most p colors for vertices of degree greater thanN/p and need not add more than ∆0 + 1 colors, where ∆0 is the maximum
degree among the remaining vertices (∆ 0 ≤ M 0 ). Overall, this means that
we have at most p + M 0 + 1 colors and hence the vector Bcolors still fits ona single processor. So, Bcolors can be initialized in a preprocessing step intime O(M 0 ).
After that, coloring any vertex v can be done in time proportional to thesize of Av , which is bounded by the degree of v. Thus, the overall time spentper processor in coloring vertices is O(M 0 ). By Lemma 5, the expected total
time (including recursions) per processor is O(M ).
Lemmas 7, 8, and 9 show that the contributions of base case, highdegree, and recurse in Algorithm 2 are within O(M ) per processor, proving
the claim on the total amount of work per processor.
As for processor idle time, observe that the bottleneck in all the algorithms as presented is the sequential processing of parts of the graphs byprocessor 1. Since the total run time (of Algorithm 3) on processor 1 is
expected to be O(M ), the same expected bound holds for the idle time ofthe other processors.

5

An add-on to achieve a good worst case behavior

So far, for a possible implementation of our algorithm, we have some degreeof freedom in choosing the number of timeslots k. If our goal is just to get
results based on expected values as shown in Section 4, we can avoid recursion
by choosing k = p and by replacing the recursive call in Algorithm 2 by acall to Sequential∆ + 1Coloring(Grec , {Fv }v ) (Algorithm 3). We can do thissince by Lemma 4 the expected size of G rec is N/k which in this case means
N/p = M , implying that Grec fits on one processor. The resulting algorithm
would have cp supersteps, for some integer c > 1; see Lemma 6.
To get a deterministic algorithm with a good worst case bound we choose
the other extreme, namely k = 2, and replace the call to the randomized
Algorithm 6, in line group vertices of Algorithm 2, by a call to Algorithms 8 and 9. This will enable us to bound the number of edges that go
into recursion. We need to distinguish between two types of edges: internaland external edges. Internal edges have both of their endpoints on the same18

Algorithm 8: Deterministically group the vertices on processor P i intok = 2 buckets.
HandleHighDegreeVertices(G,{Fv }v ,8) (see Algorithm 5);
initialize bucket[0] and bucket[1] to empty;
foreach vertex v do
determine the number of edges connecting v to bucket[0] andbucket[1], resp;
insert v in the bucket to which it has the least number of edges;
if this bucket is full thenput the remaining vertices in the other bucket;
return ;return ;processor while external edges have their endpoints on different processors.
First we argue that internal edges are handled by the call to Algorithm 8.
Then we show that external edges are handled by the subsequent call to Algorithm 9. For internal edges, the following two points need to be observed.
1. The vertices are grouped into two timeslots of about equal degree sum.
2. Most of the internal edges are good.
To achieve the first goal, Algorithm 8 first calls Algorithm 5 to get rid ofvertices with degree ≥ M/8. The constant 8 is somewhat arbitrary andcould be replaced by any constant k 0 > 2 depending on the needs of an
implementation. Algorithm 8 groups the vertices of internal edges according
to (1) and (2) above into bucket[0] and bucket[1] that will form the twotimeslots. A bucket is said to be full when the degree sum of its vertices
becomes greater than M/2.
Proposition 1 Suppose γi M of the edges on processor Pi are external (0 ≤γi M ≤ M ). Then after an application of Algorithm 8 at least ( 14 − γ2i )M of
the edges on Pi are good internal edges and each bucket has a degree sum of
at most 5M
8 .
Proof: Considering the fact that each vertex is of degree less than M/8,
the claim for the degree sum is immediate.
To see the lower bound on the number of good internal edges, consider
the bucket B that became full. The vertices in B have a degree sum of atleast M/2 and at least ( 12 − γi )M of these edges are internal. We claim thatat least half of these internal edges are good.
For the following argument, suppose that an edge is considered only
when its second endpoint is placed into a bucket. We distinguish between
19

Algorithm 9: Determine an ordering on the k = 2 buckets on eachprocessor.
foreach processor Pi do
foreach edge (v, w) do inform the processor of w about thebucket of v;
for s = 1 . . . p do0
for r, r 0 = 0, 1 do set mrr
is = 0;0foreach edge (v, w) do add 1 to mrr
is , where Ps is the processor0of w and r and r are the buckets of v and w;
0
Broadcast all values mrr
is for s = 1, . . . , p to all other processors;
inv[1] = f alse;for s = 2 to p do
A|| = 0; A = 0;
for s0 < s do
if ¬inv[s0 ] then11
A|| = A|| + m00
ss0 + mss0 ;10
A = A + m01
ss0 + mss0
else
10
A|| = A|| + m01
ss0 + mss0 ;11
A = A + m00
ss0 + mss0
if A < A|| then inv[s] = true;else inv[s] = f alse;

two types of internal edges. Early edges join vertices both of which havebeen put into a bucket before the bucket was full, while edges for which atleast one of the endpoints was placed thereafter are called late edges .
First, observe that until one of the two buckets becomes full, both buckets
have more good internal edges than bad internal edges. So, at least one halfof the early edges are good. But, notice that all the late edges that have
one endpoint in B are also good. This is the case since the second endpointof a late edge is never placed in B.
Therefore, overall, there are at least 12 ( 12 − γi )M good internal edges.
To handle the external edges we add a call to Algorithm 9 right after0
the call to Algorithm 8. This algorithm counts the number m rr
is of edges
between all possible pairs of buckets on different processors, and broadcasts
these values to all processors. Then a quick iterative algorithm is executedon each processor to ascertain as to which of the processor’s two buckets
represents the first and second timeslot.
After having decided the order in which the buckets are processed on

20

processors P1 , . . . , Pi−1 , we compute two values for processor P i : A|| thenumber of external bad edges if we would keep the numbering of the buckets
as the timeslot numbering, and A the corresponding number if we wouldinterchange them. Depending on which value is less, the order of the twobuckets of processor Pi is kept as-is or exchanged.
Using the same type of argument as the one above, we get the following
remark.
Remark 10 Algorithm 9 ensures that overall at least 21 of the external edges
are good.
Note that the above statement is true for the whole edge set, but not necessarily for the set of edges on each processor.
Proposition 2 Algorithms 8 and 9 run with linear work and communication and in a constant number of supersteps. The assignment of timeslots
they make is such that at least 14 of the edges are good.
Proof: For the number of good edges, let γ i be the fraction of externaledges of processor Pi . The total amount of good edges can now be boundedfrom below by
1
2

p
X

!

γi M +

i=1

p 
X
1
i=1

4

−

p

XM
γi
pM
N
M=
=
= .
2
4
4
4
i=1


(4)

For the complexity claim, observe that because of the load balancing done
in line balance load of Algorithm 2 (see Section 3.2.2) all processors holdthe same amount (up to a constant factor) M 0 of edges.
This implies that the recursion depth is dlog 4 pe (which is greater than
3
dlog 2 pe of the average case). Moreover, more edges go into recursion here
than in the average case and therefore the work and total communication
costs are slightly greater than the costs for the average case, but still withinO(M ).

6

Conclusion

We have presented a randomized as well as a deterministic Coarse Grained
Multicomputer coloring algorithm that colors the vertices of a general graph
G using at most ∆ + 1 colors, where ∆ is the maximum degree in G. We
showed that on a p-processor CGM model our algorithms require a paralleltime of O( |G|p ) and a total work and overall communication cost of O(|G|).
These bounds correspond to the average case for the randomized version andto the worst case for the deterministic variant. To the best of our knowledge,
our algorithms are the first parallel coloring algorithms with good speedup
for a large variety of architectures.
21

In light of the fact that LF∆+1-coloring is P-complete, a CGM LF∆+1coloring algorithm, if found, would be of significant theoretical importance.
Brent’s scheduling principle shows that anything that works well on PRAM
should, in principle, also work well on CGM or similar models (although
the constants that are introduced might be too big to be practical). But
the converse may not necessarily be true. There are P-complete problems(problems where we can’t expect exponential speedup on PRAM) that havepolynomial speedup [25]. It can be envisioned that such problems may as
well have efficient CGM algorithms. In this regard, our CGM ∆ + 1-coloringalgorithm might be a first step towards a CGM LF∆ + 1-coloring algorithm.
We also believe that, in general, designing a parallel algorithm on the
CGM model is of practical relevance. In particular, we believe that thealgorithms presented in this paper should have efficient and scalable implementations (i.e. implementations that yield good speedup for a wide rangeof N/p).
Such implementations would also be of interest in other contexts. One
example is the problem of finding maximal independent sets. Notice thatin our algorithms the vertices colored by the least color always constitute a
maximal independent set in the input graph. In general, considering G ≥i ,the subgraph induced by the color classes i, i + 1, . . ., we see that color classi always forms a maximal independent set in the graph G ≥i .

Acknowledgements
We are grateful to the anonymous referees for their helpful comments andpointing out some errors in an earlier version of this paper. Isabelle Guérin
Lassous and Jens Gustedt would like to acknowledge that part of this work
was accomplished within the framework of the “opération CGM” of theparallel research center Centre Charles Hermitte in Nancy, France.

References[1] J. R. Allwright, R. Bordawekar, P. D. Coddington, K. Dincer, and C. L.
Martin. A comparison of parallel graph coloring algorithms. Technical
Report SCCS-666, Northeast Parallel Architecture Center, Syracuse
University, 1995.
[2] G.J. Chaitin, M. Auslander, A.K. Chandra, J. Cocke, M.E. Hopkins,
and P. Markstein. Register allocation via coloring. Computer Languages, 6:47–57, 1981.
22

[3] B. Chelbus, K. Diks, W. Rytter, and T. Szymacha. Parallel complexity of lexicographically first problems for tree-structured graphs. In
A. Kreczmar and G. Mirkowska, editors, Mathematical Foundations ofComputer Science 1989, volume 379 of LNCS, pages 185–195. SpringerVerlag, 1989.
[4] T.F. Coleman and J.J. Moré. Estimation of sparse jacobian matrices
and graph coloring problems. SIAM Journal on Numerical Analysis,
20(1):187–209, 1983.
[5] D. Culler, R. Karp, D. Patterson, A. Sahay, K.E. Schauser, E. Santos,
R. Subramonian, and T. von Eicken. LogP: Towards a realistic modelof parallel computation. In Proceeding of the fourth ACM SIGPLAN
Symposium on Principles and Practice of Parallel Programming, pages
1–12, San Diego, CA, 1993.
[6] F. Dehne, A. Fabri, and A. Rau-Chaplin. Scalable parallel computational geometry for coarse grained multicomputers. International Journal on Computational Geometry, 6(3):379–400, 1996.
[7] S. Fortune and J. Wyllie. Parallelism in Random Access Machines. In
Tenth ACM Symposium on Theory of Computing, pages 114–118, San
Diego, CA, 1978.
[8] A. Gamst. Some lower bounds for a class of frequency assignment
problems. IEEE transactions of Vehicular Technology, 35(1):8–14, 1986.
[9] M.R. Garey and D.S. Johnson. Computers and Intractability. W.H.
Freeman, New York, 1979.
[10] A. H. Gebremedhin and F. Manne. Scalable parallel graph coloring algorithms. Concurrency: Practice and Experience, 12:1131–1146, 2000.
[11] M. Goudreau, K. Lang, S. Rao, T. Suel, and T. Tsantilas. Towards
efficiency and portability: Programming with the BSP model. In
8th Annual ACM symposium on Parallel Algorithms and Architectures(SPAA’96), pages 1–12, Padua, Italy, June 1996.
[12] R. Greenlaw, H.J. Hoover, and W. L. Ruzzo. Limits to Parallel Computation: P-Completeness Theory. Oxford University Press, New York,USA, 1995.
[13] I. Guérin Lassous, J. Gustedt, and M. Morvan. Feasability, portability,
predictability and efficiency: Four ambitious goals for the design and23

implementation of parallel coarse grained graph algorithms. Technical
Report 3885, INRIA, 2000.
[14] I. Guérin Lassous, J. Gustedt, and M. Morvan. Handling graphs according to a coarse grained approach: Experiments with MPI and PVM.
In Jack Dongarra, Péter Kacsuk, and N. Podhorszki, editors, Recent
Advances in Parallel Virtual Machine and Message Passing Interface,
7th European PVM/MPI Users’ Group Meeting, volume 1908 of LNCS,
pages 72–79. Springer Verlag, 2000.
[15] P. Hajnal and E. Szemerédi. Brooks coloring in parallel. SIAM Journal
on Discrete Mathematics, 3(1):74–80, 1990.
[16] M. T. Jones and P. E. Plassmann. A parallel graph coloring heuristic.
SIAM Journal of Scientific Computing, 14(3):654–669, May 1993.
[17] M. Karchmer and J. Naor. A fast parallel algorithm to color a graphwith D colors. Journal of Algorithms, 9(1):83–91, 1988.
[18] H. J. Karloff. An NC algorithm for Brook’s theorem. Theoretical Computer Science, 68(1):89–103, 16 October 1989.
[19] R. Karp and V. Ramachandran. Parallel algorithms for shared-memorymachines. In Handbook of Theoretical Computer Science Volume A:
Algorithms and Complexity, pages 869–942. Elsevier, 1990.
[20] G. Lewandowski. Practical Implementations and Applications Of Graph
Coloring. PhD thesis, University of Wisconsin-Madison, August 1994.
[21] M. Luby. A simple parallel algorithm for the maximal independent set
problem. SIAM Journal on Computing, 15(4):1036–1053, 1986.
[22] F. Manne. A parallel algorithm for computing the extremal eigenvaluesof very large sparse matrices (extended abstract). In Para98, volume1541 of LNCS, pages 332–336. Springer-Verlag, 1998.
[23] J. Naor. A fast parallel coloring of planar graphs with five colors.
Information Processing Letters, 25(1):51–53, 20 April 1987.
[24] L. G. Valiant. A bridging model for parallel computation. Communications of the ACM, 33(8):103–111, 1990.
[25] J. S. Vitter and R. A. Simons. New classes for parallel complexity: A
study of unification and other complete problems for P. IEEE Transactions on Computers, C-35(5):403–418, 1986.

24

PRO: a Model for Parallel Resource-Optimal
Computation ∗Assefaw Hadish Gebremedhin†Jens Gustedt§

Isabelle Guérin Lassous‡
Jan Arne Telle

Abstract
We present a new parallel computation model that enables the design of resource-optimal scalable parallel algorithms and simplifies theiranalysis. The model rests on the novel idea of incorporating relative
optimality as an integral part and measuring the quality of a parallel
algorithm in terms of granularity.
Key words: Parallel computers, Parallel models, Parallel algorithms, Complexity analysis∗
Research supported by IS-AUR 02-34 of The Aurora Programme, a France-Norway
Collaboration Research Project of The Research Council of Norway, The French Ministryof Foreign Affairs and The Ministry of Education, Research and Technology.
†
Department of Informatics, University of Bergen, N-5020, Norway. {assefaw,
telle}@ii.uib.no‡
LIP & INRIA Rhone-Alpes, France. Isabelle.Guerin-Lassous@inria.fr§
LORIA & INRIA Lorraine, France. gustedt@loria.fr

1

1

Introduction

One of the challenges in parallel processing is the development of a general purpose and effective model of parallel computation. Unlike the realmof sequential computation, where the Random Access Machine (RAM) has
successfully served as a standard computational model, no such single unifying model exists in the field of parallel computation. From an algorithmic point of view, the performance of a sequential algorithm is adequately
evaluated using its execution time making the RAM powerful enough foranalysis and design. On the other hand, the performance evaluation of aparallel algorithm involves several metrics, the most important of which are
speedup, optimality (or efficiency), and scalability. Speedup and optimality
are relative in nature as they are expressed with respect to some sequential
algorithm. The notion of relativity is also relevant from a practical pointof view. A parallel algorithm is often not designed from scratch, but rather
starting from a sequential algorithm.
We believe that a parallel computation model should incorporate themost important performance evaluation metrics of parallel algorithms as
the RAM does for sequential algorithms. In light of this, the objective of
the current work is to develop a model that simplifies the design and analysisof resource-optimal scalable parallel algorithms.
In an interesting survey paper [21], Maggs et al. suggest that an ideal
parallel computation model be designed within “the philosophy of simplicity and descriptivity balanced with prescriptivity”. The Parallel ResourceOptimal (PRO) computation model proposed here is developed within thisspirit. The key features of the PRO model that distinguish it from existing
parallel computation models are relativity, resource-optimality, and a new
quality measure referred to as granularity.
Relativity pertains to the fact that the design and analysis of a parallel
algorithm in PRO is done relative to the time and space complexity of a
specific sequential algorithm. Consequently, the parameters involved in theanalysis of a PRO-algorithm are the number of processors p, the input sizen, and the time and space complexity of the reference sequential algorithm
Aseq .
A PRO-algorithm is required to be both time- and space-optimal (hence
resource-optimal). A parallel algorithm is said to be time- (or work-) optimal
if the overall computation and communication cost involved in the algorithm
is proportional to the time complexity of the sequential algorithm used asa reference. Similarly, it is said to be space-optimal if the overall memory
space used by the algorithm is of the same order as the memory usage of2the underlying sequential version. As a consequence of its time-optimality, a
PRO-algorithm always yields linear speedup relative to the reference sequential algorithm; i.e., the ratio between the sequential and parallel runtime isa linear function of p.
The quality of a PRO-algorithm is measured by the range of values p
can assume while linear speedup is maintained. This range is capturedby an attribute of the model called the granularity function Grain(n). In
other words, a PRO-algorithm with granularity Grain(n) is required to be
fully scalable for all values of p such that p = O(Grain(n)). The granularity
function Grain(n) determines the quality of one PRO-algorithm over anotherrelative to the same sequential time and space complexity. The higher thefunction value Grain(n) the better the algorithm. Note that since optimality(consequently linear speedup) is ‘hard-wired’ into the model, the runtime
cannot be a quality measure for a PRO algorithm. However, in a sense, thetime and space complexity of the reference sequential algorithm A seq can also
be seen as a quality measure of the PRO-algorithm. This means that the
selection of the reference sequential algorithm is of significant importance.
The rest of the paper is organized as follows. In Section 2 we give an
overview of existing parallel computation models and highlight their limitations. In Section 3 the PRO model is presented in detail and in Section 4
it is compared with a selection of existing parallel models. In Section 5 we
illustrate how the model is used in design and analysis using the matrix multiplication problem as an example. In Section 6 we give a PRO-algorithm
for one-to-all broadcast, as an example of a primitive communication routine found in a potential PRO library. Finally, we conclude the paper in
Section 7 with some remarks.

2

Existing models and their limitations

There exists a plethora of parallel computation models in the literature. On
the theoretical end, we find the Parallel Random Access Machine (PRAM)
model [8, 17] which in its simplest form posits a set of p processors, with
global shared memory, executing the same program in lockstep. In thismodel, every processor can access any memory location at unit cost of time
regardless of the memory location. This assumption is in obvious disagreement with the reality of practical parallel computers.
However, despite its serious limitation of being an ‘idealized’ model ofparallel computation, the standard PRAM model still serves as a theoretical framework for investigating the maximum possible computational parallelism in a given task. Specifically, on this model, the N C versus P -complete3dichotomy [14] is used to reflect the ease/hardness of finding a parallel algorithm for a problem. Recall that N C denotes the class of problems whichhave PRAM-algorithms with polylogarithmic runtime and polynomial number of processors in the input size. A problem is said to be P -complete
if an N C-algorithm for it would imply that all polynomial time sequential
problems have N C-algorithms. The problem of whether or not P = N C
has long been an open problem.
The N C versus P -complete dichotomy has its own practical limitations.
First, P -completeness does not depict a full picture of non-parallelizability
since the runtime requirement for an N C parallel algorithm is so stringent
that the classification is confined to the case where up to polynomial number of processors in the input size is available (fine-grained setting). For
example, there are P -complete problems for which less ambitious, but stillsatisfactory, runtime can be obtained by parallelization in PRAM [23]. Ina fine-grained setting, since the number of processors p is a function of the
input size n, it is customary to express speedup as a function of n. Thus
the speedup obtained using an N C-algorithm is sometimes referred to as
exponential. In a coarse-grained setting, i.e., the case where n and p areorders of magnitude apart, speedup is expressed as a function of only p andsome recent results [4, 7, 9, 15] show that this approach is practically relevant. Second, an N C-algorithm is not necessarily work-optimal, and thus
not resource-optimal considering runtime and memory space as resourcesthat one wants to use efficiently. Third, even if we restrict ourselves towork-optimal N C-algorithms and apply Brent’s scheduling principle, whichsays an algorithm in theory can be simulated on a machine with fewer processors by only a constant factor more work, implementations of PRAM
algorithms often do not reflect this optimality in practice [6]. This is mainly
because the PRAM model does not account for non-local memory access(communication), and a Brent-type simulation relies heavily on cheap communication.
To overcome the defects of the PRAM related to its failure of capturing real machine characteristics, the advocates of shared memory models
propose several modifications to the standard PRAM model. In particular,
they enhance the standard PRAM model by taking practical machine features such as memory access, synchronization, latency and bandwidth issues
into account. Pointers to the PRAM family of models can be found in [21].
Critics of shared memory models argue that the PRAM family of models fail to capture the nature of existing parallel computers with distributed
memory architectures. Examples of distributed memory computational models suggested as alternatives include the Postal Model [2] and the Block
4

Distributed Memory (BDM) model [18]. Other categories of parallel models such as low-level, hierarchical memory, and network models are briefly
reviewed in [21].
A more recent category of parallel models is that of ‘bridging’ models, a
notion popularized by Valiant with his introduction of the Bulk Synchronous
Parallel (BSP) model [22]. The BSP model is a distributed memory coarsegrained model in which parallel computation proceeds as a sequence of barrier synchronized supersteps where local computation and communication
are distinct rather than intermingled phases. Culler et al. [5] extended the
BSP model by allowing asynchronous execution and better accounting forcommunication overhead. Their model is coined LogP, an acronym for thefour parameters involved. A common feature of the BSP, LogP, and otherrelated models is their lack of simplicity: each model involves relatively
many parameters making analysis and design of algorithms cumbersome.
The Coarse Grained Multicomputer (CGM) model [4, 7] was later proposed in an effort to retain the advantages of BSP while keeping the model
simple (making the number of parameters fewer). The BSP and its special
case CGM have been the primary inspirations for our model. Thus, we believe that many optimal CGM and BSP algorithms can easily be adapted
to PRO.
The PRO model attempts to partially address the limitations of existing
parallel models highlighted in the foregoing discussion and compromises between theoretical and practical considerations. One of its advantages froma theoretical point of view is that it is a step forward towards the identification of the class of problems for which ‘good’ parallel algorithms exist
in a more realistic (practical) way than the existing N C versus P -complete
classification.
Our main goal in suggesting the PRO model is to enable the developmentof scalable and resource-optimal parallel algorithms and to simplify theiranalysis. The model identifies the salient features of a parallel algorithm thatmake its practical scalability and optimality highly likely. In this regard, it
can be considered as a set of ‘guidelines’ for the algorithm designer in thequest for developing scalable and efficient parallel algorithms. Hence, PRO
can be seen as a mix of a parallel computation model and a parallel algorithm
design scheme which makes it biased towards the software side in its role asa bridging model.

5

3

The PRO model

The PRO model is an algorithm design and analysis tool used to delivera practical, optimal, and scalable parallel algorithm relative to a specificsequential algorithm whenever this is possible. Let Time(n) and Space(n)denote the time and space complexity of a specific sequential algorithm fora given problem with input size n. The PRO model is defined to have thefollowing attributes.
Machine The underlying machine is assumed to consist of p processors
with M = O( Space(n)) private memory each, interconnected by somepcommunication network (or shared memory) that can deliver messages
in a point-to-point fashion. A message can consist of several machinewords.
Coarseness We assume that p ≤ M , i.e., the size of the local memory ofeach processor is big enough to store p words.
Execution For any value p = O(Grain(n)), a PRO algorithm,• consists of O( Time(n)) supersteps. A superstep consists of a local
p2
computation phase and an interprocessor communication phase.
In particular, in each superstep, each processor– sends at most one message to every other processor,– sends and receives at most M words in total, and pays a unitof time per word sent and received,– performs local computation, and pays a unit of time per operation,).
• has parallel runtime Time(n, p) = O( Time(n)pNote that the granularity function Grain(n) is a quality measure of a
PRO-algorithm.
As discussed in the LogP paper [5], technological factors are forcing parallel systems to converge towards systems formed by a collection of essentially complete computers connected by a robust communication network.
The machine model assumption of PRO is consistent with this convergenceand maps well on several existing parallel computer architectures. The mem) ensures that the space utilized by the unory requirement M = O( Space(n)pderlying sequential algorithm is uniformly distributed among the p processors. Since we may, without loss of generality, assume that Space(n) = Ω(n),the implication is that the private memory of each processor is large enough6

to store its ‘share’ of the input and any additional space the sequential algorithm might require. When Space(n) = Θ(n), note that the input data
must be uniformly distributed on the p processors. In this case the machinemodel assumption of PRO is similar to the assumption in the CGM model[7].
The coarseness assumption p ≤ M is consistent with the structure ofexisting parallel machines and machines to be built in the foreseeable future.
The assumption is required to simplify the implementation of collecting
messages (from possibly all other processors) on a single processor.
The execution of a PRO-algorithm consists of a sequence of supersteps(or rounds). The length of (time spent in) a superstep on each processor isdetermined by the sum of the time used for communication and the timeused for local computation. The length of a superstep s in the parallel
algorithm seen as a whole, denoted by Time s (n, p), is the maximum overthe lengths of the superstep on all processors. We can conceptually think as
if the supersteps are synchronized by a barrier set at the end of the longest
superstep across the processors. However, note that in PRO the processors
are not in reality required to synchronize at the end of each superstep. The
parallel runtime Time(n, p) of the algorithm is the sum of the lengths of
all the supersteps. Notice that the hypothetical barriers result in only aconstant factor more time compared with an analysis that does not assume
the barriers.
In PRO, since a processor sends at most one message to every otherprocessor in each superstep, each processor is involved in at most 2(p − 1)messages per superstep. Therefore, the requirement Steps = O( Time(n)) onp2
the number of supersteps implies that the overall time paid per processor
for communication overhead and latency is O(Time(n)/p) and hence can beneglected from the analysis since our goal is to achieve an O(Time(n)/p)parallel runtime. Notice that the bandwidth restriction of the underlyingarchitecture which in turn contributes to the communication cost is accounted for since each processor pays a unit of time per word sent andreceived. This is not an unrealistic assumption noting that the network
throughput (accounted in machine words) on modern architectures such as
high performance clusters is relatively close to the CPU frequency and to
the CPU/memory bandwidth.
The condition Time(n, p) = O( Time(n)) requires that a PRO-algorithmpbe optimal and yield linear speedup relative to the sequential algorithm used
as a reference. This requirement ensures the potential practical use of theparallel algorithm.

7

Observation 1 A PRO algorithm relative to a sequential algorithm withruntime O(Time(n)) and space requirement O(Space(n)) has maximum granpppularity Grain(n) = O(min{ Space(n), (Time(n)}) = O( Space(n)). A
PRO algorithm that achieves this is said to have optimal grain.
Observation 1 is due to the limit on the memory size of each processor,
the coarseness assumption, and the bound on the number of supersteps. The))
limit on the size of the private memory of each processor (M = O( Space(n)pptogether with the coarseness assumption p ≤ M imply p = O( Space(n)).
The fact that the number of supersteps of a PRO-algorithm should bep
Steps = O(Time(n)/p2 ), gives p = O( (Time(n)/Steps)) upon resolving
and we clearly have Steps ≥ 1. Finally, note that Time(n) ≥ Space(n),since an algorithm has to at least read the input.
Since a PRO-algorithm yields linear speedup for any p = O(Grain(n)),
a result like Brent’s scheduling principle is implicit for these values of p.
But Observation 1 shows that we cannot start with an arbitrary number ofprocessors and efficiently simulate on a fewer number. So Brent’s schedulingprinciple does not hold with full generality in the PRO model, which is inaccordance with practical observations.
The design of a PRO-algorithm may sometimes involve subroutines forwhich there do not exist sequential counterparts. Examples of such tasks include communication primitives such as broadcasting, data (re)-distribution
routines, and load balancing routines. Such routines are often required invarious parallel algorithms. With a slight abuse of notation, we call suchparallel routines PRO-algorithms if the overall computation and communication cost is linear in the input size to the routines.

4Comparison with other models

In this section we compare the PRO model with PRAM, QSM, BSP, LogP,
and CGM. Our tabular format for comparison is inspired by a similar presentation in [13], where the Queuing Shared Memory (QSM) model is proposed.
The columns of Table 1 are labeled with the names of the selected models in
our comparison and some relevant features of a model are listed along therows.
The synchrony assumption of the model is indicated in the row labeled
synch. Lock-step indicates that the processors are fully synchronized at eachstep (of a universal clock), without accounting for synchronization. Bulksynchrony indicates that there can be asynchronous operations between synchronization barriers. The row labeled memory shows how the model views
8synch.
memorycommun.
parameters
granularity
speedup
optimalquality

PRAM [8]lock-step
sh.
SM
n
fine
NA
NA
time

QSM [13]bulk-synch.
sh.
SMp, g, n
fine
NA
NA
time

BSP [22]bulk-synch.
dist.
MPp, g, L, ncoarseNA
NA
time

LogP [5]asynch.
dist.
MPp, g, l, o, n
fine
NA
NA
time

CGM [4]
asynch.
priv.
MP/SM
p, n
coarseNANAroundsPRO
asynch.
priv.
MP/SMp, n, Aseq
Grain(n)Θ(p)
rel. Aseq
Grain(n)Table 1: Comparison of parallel computational models
the memory of the parallel computer: sh. indicates globally accessible sharedmemory, dist. stands for distributed memory and priv. is an abstraction
for the case where the only assumption is that each processor has accessto private (local) memory. In the last variant the whole memory could either be distributed or shared. The row labeled commun. shows the typeof interprocessor communication assumed by the model. Shared memory(SM) indicates that communication is effected by reading from and writing
to a globally accessible shared memory. Message-passing (MP) denotes thesituation where processors communicate by explicitly exchanging messages
in a point-to-point fashion. The MP abstraction hides the details of how
the message is routed through the interprocessor communication network.
The parameters involved in the model are indicated in the row labeledparameters. The number of processors is denoted by p, n is the input size,
Aseq is the reference sequential algorithm, l is the communication cost (latency), L is a single parameter that accounts for the sum of latency (l)and the cost for a barrier synchronization, g is the bandwidth gap, and ois the overhead associated with sending or receiving a message. Note thatthe machine characteristics l and o are are taken into account in PRO, eventhough they are not explicitly used as parameters. Latency is taken into
consideration since the length of a superstep is determined by the sum of
the computational and communication cost. Communication overhead is
hidden by the PRO-requirement that states Steps = O( Time(n)).
p2
The row labeled granularity indicates whether the model is fine-grained,coarse-grained or a more precise measure is used. We say that a model iscoarse-grained if it applies to the case where n  p and call it fine-grained if
it relies on using up to a polynomial number of processors in the input size.
In PRO granularity is exactly the quality measure Grain(n), and appears as
one of the attributes of the model.
The rows labeled speedup and optimal indicate the speedup and resource
optimality requirements imposed by the model. Whenever these issues are
not directly addressed by the model or are not applicable, the word ‘NA’ is9used. Note that these requirements are ‘hard-wired’ in the model in the caseof PRO. The label ‘rel. Aseq ’ means that the algorithm is optimal relative
to the time and space complexity of A seq . We point out that the goal in thedesign of algorithms using the CGM model [7, 4] is usually stated as thatof achieving optimal algorithms, but the model per se does not impose an
optimality requirement.
The last row indicates the quality measure of an algorithm designed
using the different models. For all other models except CGM and PRO,
the quality measure is running time. In CGM, the number of supersteps(rounds) is usually presented as a quality measure. In PRO the quality
measure is granularity, one of the features that make PRO fundamentallydifferent from all existing parallel computation models.

5

Algorithm example: matrix multiplication

In this section we illustrate how the PRO model is used, by starting from agiven sequential algorithm and then designing and analyzing a parallel algorithm relative to it. We use the standard matrix multiplication algorithmwith three nested for-loops as an example. This example is chosen for its
simplicity and since our objective at this stage is to illustrate the use of a
new model rather than solving a “difficult” problem.
Consider the problem of computing the product C of two m × m matrices A and B (input size n = m2 ). We want to design a PRO-algorithm
relative to the standard sequential matrix multiplication algorithm which
3
has Time(n) = O(n 2 ) and Space(n) = O(n).
We assume that the input matrices A and B are distributed among
the p processors P0 , ..., Pp−1 so that processor Pi stores rows (respectivelymcolumns) m
p · i + 1 to p · (i + 1) of A (respectively B). The output matrix C
will be row-partitioned among the p processors in a similar fashion. Noticethat with this data distribution each processor can, without communica2
m2
tion, compute a block of mp2 of the p entries of C expected to reside on
2

it. In order to compute the next block of mp2 entries, processor Pi needs
the columns of matrix B that reside on processor P i+1 . In each superstep
the processors in the PRO algorithm will therefore exchange columns in around-robin fashion and then each will compute a new block of results. Notethat each column exchanged in a superstep constitutes one single message.
Note also that the initial distribution of the rows of matrix A remains unchanged. In Algorithm 1, we have organized this sequence of computationand communication steps in a manner that meets the requirements of the10

Algorithm 1: Matrix multiplication
Input: Two m × m matrices A and B. The rows (columns) of A
(B) are divided into m/p contiguous blocks, and stored onprocessors P0 , P1 , . . . Pp−1 respectively
Output: The product matrix C where the rows are stored in contiguous blocks across the p processors
for superstep s = 1 to p do
foreach processor Pi do
Pi computes the local sub-matrix product of its rows andcurrent columns;
P(i+1)mod p sends its current block of columns to P i ;Pi receives a new current block of columns from P (i+1)mod p ;PRO model.
Algorithm 1 has p supersteps (Steps = p). In each superstep, the timespent in locally computing each of the m 2 /p2 entries is Θ(m) resulting in
3
local computing time Θ(m3 /p2 ) = Θ(n 2 /p2 ) per superstep. Likewise, thetotal size of data (words) exchanged by each processor in a superstep isΘ(m2 /p) = Θ(n/p). Thus, the length of a superstep s is Time s (n, p) =√3
3
Θ(n 2 /p2 +n/p). Note that for p = O( n), Times (n, p) = Θ(n 2 /p2 ). Hence,√for p = O( n), the overall parallel runtime of the algorithm isTime(n, p) =

X

3

3

Θ(n 2 /p2 ) = Θ(n 2 /p) = Θ(Time(n)/p).

(1)

Steps

Noting that Space(n) = Θ(n), we see that the memory restriction of
the PRO model is respected, i.e., each processor has enough memory size
to handle the transactions. In order to be able to neglect communication
overhead, the condition on the number of supersteps, which in this caseis just p, should be met. In other words, we need p = O(Time(n)/p 2 ) =√3
O(n 2 /p2 ), which is true for p = O( n). Thus the granularity function of√
the PRO-algorithm is Grain(n) = n.
In summary,
Lemma 1 Multiplication of two m by m matrices has a PRO-algorithmwith Grain(n) = m relative to a sequential algorithm with Time(n) = m 3
and Space(n) = m2 (input size n = m2 ).
From Observation 1, we note that Algorithm 1 achieves optimal granularity. Note that on a relaxed model, where the assumption that p ≤ M
is not present, the strong regularity of matrix multiplication and the exact
11

knowledge of the communication pattern allows for algorithms that have aneven finer granularity than m. For example, a systolic matrix multiplication algorithm has a granularity of m 2 . However, PRO is intended to be
applicable for general problems and practically relevant parallel systems.

6

Communication primitive example: one-to-all
broadcast

A good parallel computation model should have a selection of algorithms
for primitive communication tasks available in its algorithm design toolbox. The PRO model is intended to meet this demand, but for lack of space
we give only one example.
In this section we illustrate how the PRO model allows optimal oneto-all broadcasting among its processors. Since there is no sequential basis
algorithm in this case, we want an algorithm whose overall communicationand computation cost is linear in the input and output sizes. More precisely,
we consider the situation where the input consists of a vector of size m ona single processor and the output should be a copy of this vector on each of
the p processors, and we want an algorithm that achieves this in O(m) timeusing O(m) memory on each processor. See Algorithm 2.
Algorithm 2: One-to-All Broadcast
Input: A vector V of size m on processor P 0
Output: A copy of V on each processor
S1 P0 divides V into p equal sized parts;
P0 sends the ith part of V to processor Pi , for each 0 < i ≤ p;foreach processor Pi , i > 0 do processor Pi receives the ith partfrom P0 ;S2 foreach processor Pi do
Pi sends out the ith part to Pj , for each j 6= i and 0 < j ≤ p.
foreach processor Pj ,j 6= 0 do
Pj receives the ith part from Pi , for each i 6= j and 0 < i ≤ p

Lemma 2 PRO Algorithm 2 implements a one-to-all broadcast of m memory words in two supersteps using O(m) time and O(m) space per processor,
for any number of processors p ≤ m.
Proof: First, we note that the algorithm correctly broadcasts the desired
vector V , while observing the space restriction, in two supersteps. We turnto the timing. In step S1 processor P 0 in total sends out (p − 1)m/p words
12and each of the other processors receives a message of size m/p. In step S2
processor Pi in total sends out p−2
p m words. Processor Pj , j 6= 0, in total
p−1
receives p m words.
The total time is dominated by the communication which is
p−2
p−1
m+
m=
p
p
m/p(p + p − 2 + p − 1) < 3m(p − 1)m/p + m/p +for total time O(m) as claimed.

7

(2)
(3)
2

Conclusion

We have introduced a new parallel computation model (called PRO) that enables the development of efficient scalable parallel algorithms and simplifies
the complexity analysis of such algorithms.
The distinguishing feature of the PRO model is the novel focus on relativity, resource-optimality, and a new quality measure (granularity). In
particular, the model requires a parallel algorithm to be both time- andspace-optimal relative to an underlying sequential algorithm. Having optimality as a built-in requirement, the quality of a PRO-algorithm is measured
by the maximum number of processors that could be used while the optimality of the algorithm is maintained.
The focus on relativity has theoretical as well as practical justifications.
From a theoretical point of view, the performance evaluation metrics ofa parallel algorithm includes speedup and optimality, both of which arealways expressed relative to some sequential algorithm. Moreover, there
is an inherent asymmetry between sequential and parallel computation. A
parallel algorithm would always imply a sequential algorithm, whereas theconverse is usually not true. Thus, in a sense, it is natural to think of anunderlying sequential algorithm whenever one speaks of a parallel algorithm.
From a practical point of view, one notes that the development of a parallel
algorithm is often built on some known sequential algorithm.
The fact that optimality is incorporated as a requirement in the PROmodel enables one to concentrate only on parallel algorithms that are practically useful.
However, the PRO model is not just a collection of some ‘ideal’ features ofparallel algorithms, it is also a means to achieve these features. In particular,
the attributes of the model capture the salient characteristics of a parallel
algorithm that make its practical optimality and scalability highly likely.
13

In this sense, it can also be seen as a parallel algorithm design scheme.
Moreover, the simplicity of the model eases analysis.
We believe that the PRO model is a step forward towards the identification of problems for which ‘practically good’ parallel algorithms exist.
Much work remains to be done, and we hope that other members of theresearch community will join in. As a first item on the agenda, the PROmodel needs to be tested for compatibility with already existing practicalparallel algorithms.
Acknowledgments
helpful comments.

We are grateful to the anonymous referees for their

References[1] A. G. Alexandrakis, A. V. Gerbessiotis, D. S. Lecomber, and C. J.
Siniolakis. Bandwidth, space and computation efficient PRAM programming: The BSP approach. In Proceedings of the SUP’EUR ’96
Conference, Krakow, Poland, September 1996.
[2] A. Bar-Noy and S. Kipnis. Designing broadcasting algorithms in the
Postal Model for message passing systems. In The 4th annual ACM
symposium on parallel algorithms and architectures, pages 13–22, July
1992.
[3] R. P. Brent. The parallel evaluation of generic arithmetic expressions.
Journal of the ACM, 21(2):201–206, 1974.
[4] E. Caceres, F. Dehne, A. Ferreira, P. Locchini, I. Rieping, A. Roncato,
N. Santoro, and S. W. Song. Efficient parallel graph algorithms forcoarse grained multicomputers and BSP. In The 24th International
Colloquium on Automata Languages and Programming, volume 1256 of
LNCS, pages 390–400. Springer Verlag, 1997.
[5] D. E. Culler, R. M. Karp, D. A. Patterson, A. Sahay, K. E. Schauser,
E. Santos, R. Subramonian, and T. von Eicken. LogP: Towards a realistic model of parallel computation. In 4th ACM SIGPLAN Symposium
on principles and practice of parallel programming, San Diego, CA, May
1993.
[6] F. Dehne. Coarse grained parallel algorithms. Algorithmica Special
Issue on “Coarse grained parallel algorithms”, 24(3/4):173–176, 1999.

14

[7] F. Dehne, A. Fabri, and A. Rau-Chaplin. Scalable parallel computational geometry for coarse grained multicomputers. International Journal on Computational Geometry, 6(3):379–400, 1996.
[8] S. Fortune and J. Wyllie. Parallelism in random access machines. In
10th ACM Symposium on Theory of Computing, pages 114–118, May
1978.
[9] A. H. Gebremedhin, I. Guérin Lassous, J. Gustedt, and J. A. Telle.
Graph coloring on a coarse grained multiprocessor. In Ulrik Brandes
and Dorothea Wagner, editors, WG 2000, volume 1928 of LNCS, pages
184–195. Springer-Verlag, 2000.
[10] A. V. Gerbessiotis, D. S. Lecomber, C. J. Siniolakis, and K. R. Sujithan. PRAM programming: Theory vs. practice. In Proceedings of 6thEuromicro Workshop on Parallel and Distributed Processing, Madrid,
Spain. IEEE Computer Society Press, January 1998.
[11] A. V. Gerbessiotis and C. J. Siniolakis. A new randomized sorting
algorithm on the BSP model. Technical report, New Jersey Institute ofTechnology, 2001.
[12] A. V. Gerbessiotis and L. G. Valiant. Direct bulk-synchronous parallel
algorithms. Journal of Parallel and Distributed Computing, 22:251–267,
1994.
[13] P. B. Gibbons, Y. Matias, and V. Ramachandran. Can a SharedMemory Model Serve as a Bridging Model for Parallel Computation?
Theory of Computing Systems, 32(3):327–359, 1999.
[14] R. Greenlaw, H.J. Hoover, and W. L. Ruzzo. Limits to Parallel Computation: P-Completeness Theory. Oxford University Press, New York,
1995.
[15] I. Guérin Lassous, J. Gustedt, and M. Morvan. Handling graphs according to a coarse grained approach: Experiments with MPI and PVM.
In Jack Dongarra, Péter Kacsuk, and N. Podhorszki, editors, 7th European PVM/MPI Users’ Group Meeting, volume 1908 of LNCS, pages
72–79. Springer Verlag, 2000.
[16] K. Hawick et al. High performance computing and communications
glossary. see http://nhse.npac.syr.edu/hpccgloss/.
[17] J. Jájá. An Introduction to Parallel Algorithms. Addison-Wesley, 1992.
15

[18] J. JáJá and K. W. Ryu. The Block Distributed Memory model. IEEE
Transactions on Parallel and Distributed Systems, 8(7):830–840, 1996.
[19] R. M. Karp and V. Ramachandran. Parallel Algorithms for SharedMemory Machines. In Jan van Leeuwen, editor, Handbook of TheoreticalComputer Science, volume A, Algorithms and Complexity, pages 869–
941. Elsevier Science Publishers B.V., Amsterdam, 1990.
[20] C. P. Kruskal, L. Rudolph, and M. Snir. A complexity theory of efficientparallel algorithms. Theoretical Computer Science, 71(1):95–132, march
1990.
[21] B. M. Maggs, L. R. Matheson, and R. E. Tarjan. Models of parallel
computation: A survey and synthesis. In 28th HICSS, volume 2, pages
61–70, January 1995.
[22] L. G. Valiant. A bridging model for parallel computation. Communications of the ACM, 33(8):103–111, 1990.
[23] J. S. Vitter and R. A. Simons. New classes for parallel complexity: A
study of unification and other complete problems for P. IEEE Transactions on Computers, C-35(5):403–418, 1986.

16