c 2005 Society for Industrial and Applied Mathematics


SIAM REVIEW
Vol. 47, No. 4, pp. 629–705

What Color Is Your Jacobian?
Graph Coloring for
Computing Derivatives∗
Assefaw Hadish Gebremedhin†
Fredrik Manne‡Alex Pothen†Abstract. Graph coloring has been employed since the 1980s to eﬃciently compute sparse Jacobianand Hessian matrices using either ﬁnite diﬀerences or automatic diﬀerentiation. Several
coloring problems occur in this context, depending on whether the matrix is a Jacobian or a
Hessian, and on the speciﬁcs of the computational techniques employed. We consider eightvariant vertex coloring problems here. This article begins with a gentle introduction to the
problem of computing a sparse Jacobian, followed by an overview of the historical development of the research area. Then we present a unifying framework for the graph modelsof the variant matrix estimation problems. The framework is based upon the viewpointthat a partition of a matrix into structurally orthogonal groups of columns corresponds
to distance-2 coloring an appropriate graph representation. The uniﬁed framework helps
integrate earlier work and leads to fresh insights; enables the design of more eﬃcient algorithms for many problems; leads to new algorithms for others; and eases the task ofbuilding graph models for new problems. We report computational results on two of thecoloring problems to support our claims. Most of the methods for these problems treat acolumn or a row of a matrix as an atomic entity, and partition the columns or rows (orboth). A brief review of methods that do not ﬁt these criteria is provided. We also discuss
results in discrete mathematics and theoretical computer science that intersect with the
topics considered here.
Key words. sparsity, symmetry, Jacobians, Hessians, ﬁnite diﬀerences, automatic diﬀerentiation, matrix partitioning problems, distance-k coloring, approximation algorithms
AMS subject classiﬁcations. 05C15, 05C90, 90C06, 90C27, 90C90
DOI. 10.1137/S0036144504444711

When colour goes home into the eyes,
And lights that shine are shut again . . .
And that no-place which gave them birth, shall close
The rainbow and the rose . . .
—Rupert Brooke, “The Treasure”∗ Received by the editors July 1, 2004; accepted for publication (in revised form) April 4, 2005;
published electronically October 31, 2005.
http://www.siam.org/journals/sirev/47-4/44471.html† Computer Science Department, Old Dominion University, Norfolk, VA 23529-0162 (assefaw@
cs.odu.edu, pothen@cs.odu.edu). The work of the ﬁrst author was supported by the Universityof Bergen while completing his Ph.D. and later by U.S. National Science Foundation grant ACI
0203722. The work of the third author was supported by U.S. National Science Foundation grant
ACI 0203722, by DOE grant DE-FC02-01ER25476, by subcontracts B533520 and B542604 from theLawrence Livermore National Laboratory, and by the Computer Science Research Institute of theSandia National Laboratories during a sabbatical year.
‡ Department of Informatics, University of Bergen, N-5020 Bergen, Norway (fredrikm@ii.uib.no).

629

630

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

Contents.
1 Introduction
632
1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 632
1.2 A Jacobian Computation Problem . . . . . . . . . . . . . . . . . . . . 632
1.3 Automatic Diﬀerentiation . . . . . . . . . . . . . . . . . . . . . . . . . 635
1.4 Variations on Matrix Computation . . . . . . . . . . . . . . . . . . . . 636
1.4.1 Jacobian vs. Hessian . . . . . . . . . . . . . . . . . . . . . . . . 637
1.4.2 Direct vs. Substitution-Based Evaluation . . . . . . . . . . . . 637
1.4.3 Unidirectional vs. Bidirectional Partition . . . . . . . . . . . . 638
1.4.4 Full vs. Partial Computation . . . . . . . . . . . . . . . . . . . 639
1.4.5 Other Variations . . . . . . . . . . . . . . . . . . . . . . . . . . 639
1.5 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 640
1.6 New Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 641
1.7 Scope and Further Reading . . . . . . . . . . . . . . . . . . . . . . . . 641
1.8 Terminology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 642
2 Deﬁnitions of Graph-Theoretic Concepts
642
2.1 Preliminary Concepts and Notations . . . . . . . . . . . . . . . . . . . 643
2.2 Distance-k Graph Coloring . . . . . . . . . . . . . . . . . . . . . . . . 643
2.3 The Power of a Graph . . . . . . . . . . . . . . . . . . . . . . . . . . . 644
2.4 Representing Matrices Using Graphs . . . . . . . . . . . . . . . . . . . 644
3 Unidirectional, Direct Computation of the Jacobian
645
3.1 The Matrix Partitioning Problem . . . . . . . . . . . . . . . . . . . . . 646
3.2 A Graph Coloring Formulation . . . . . . . . . . . . . . . . . . . . . . 646
3.3 An Alternative Coloring Formulation . . . . . . . . . . . . . . . . . . . 647
3.4 Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 648
3.4.1 Flexibility . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 648
3.4.2 Uniﬁcation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 648
3.4.3 Graph Size . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 648
3.4.4 Ease of Construction . . . . . . . . . . . . . . . . . . . . . . . . 649
3.4.5 Use of Existing Software . . . . . . . . . . . . . . . . . . . . . . 649
3.5 Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 649
3.5.1 A Distance-2 Coloring Algorithm . . . . . . . . . . . . . . . . . 650
3.5.2 A Partial Distance-2 Coloring Algorithm . . . . . . . . . . . . . 651
3.5.3 A Distance-1 Coloring Algorithm . . . . . . . . . . . . . . . . . 652
3.6 Experimental Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 653
3.6.1 Test Matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . 653
3.6.2 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . 653
4 Direct Computation of the Hessian658
4.1 The Matrix Partitioning Problem . . . . . . . . . . . . . . . . . . . . . 658
4.2 A Graph Coloring Formulation . . . . . . . . . . . . . . . . . . . . . . 659
4.3 Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 661
4.3.1 The First Star Coloring Algorithm . . . . . . . . . . . . . . . . 661
4.3.2 The Second Star Coloring Algorithm . . . . . . . . . . . . . . . 662
4.4 Experimental Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 663
4.4.1 Test Graphs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 664
4.4.2 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . 664

GRAPH COLORING FOR COMPUTING DERIVATIVES

631

5 Bidirectional, Direct Computation of the Jacobian664
5.1 The Need for a Bidirectional Partition . . . . . . . . . . . . . . . . . . 665
5.2 The Matrix Partitioning Problem . . . . . . . . . . . . . . . . . . . . . 665
5.3 A Graph Coloring Formulation . . . . . . . . . . . . . . . . . . . . . . 667
5.4 Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 668
6 Substitution Methods
670
6.1 Computing the Hessian . . . . . . . . . . . . . . . . . . . . . . . . . . 670
6.1.1 The Matrix Partitioning Problem . . . . . . . . . . . . . . . . . 671
6.1.2 A Graph Coloring Formulation . . . . . . . . . . . . . . . . . . 671
6.1.3 An Acyclic Coloring Algorithm . . . . . . . . . . . . . . . . . . 674
6.2 Computing the Jacobian . . . . . . . . . . . . . . . . . . . . . . . . . . 676
6.2.1 The Matrix Partitioning Problem . . . . . . . . . . . . . . . . . 676
6.2.2 A Graph Coloring Formulation . . . . . . . . . . . . . . . . . . 676
6.2.3 An Acyclic Bicoloring Algorithm . . . . . . . . . . . . . . . . . 676
7 Interrelationships among the Coloring Problems
677
7.1 Chromatic Numbers . . . . . . . . . . . . . . . . . . . . . . . . . . . . 677
7.2 Two-Colored Induced Subgraphs . . . . . . . . . . . . . . . . . . . . . 677
8 Partial Matrix Computation
679
8.1 Unidirectional Computation of the Jacobian . . . . . . . . . . . . . . . 679
8.2 Computing the Hessian . . . . . . . . . . . . . . . . . . . . . . . . . . 681
8.3 Bidirectional Computation of the Jacobian . . . . . . . . . . . . . . . . 682
9 Hypergraph Coloring Formulations
684
9.1 Deﬁnitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 684
9.2 Hypergraph Coloring . . . . . . . . . . . . . . . . . . . . . . . . . . . . 684
9.3 Representing Matrices Using Hypergraphs . . . . . . . . . . . . . . . . 685
9.4 Structurally Orthogonal Partition and Hypergraph Coloring . . . . . . 685
10 Other Matrix Estimation Methods
686
10.1 Methods Based on Solving a Rectangular System . . . . . . . . . . . . 687
10.2 Element and Variable Isolation . . . . . . . . . . . . . . . . . . . . . . 687
10.3 Element Isolation and Edge Coloring . . . . . . . . . . . . . . . . . . . 688
10.4 Partitioning Segmented Columns . . . . . . . . . . . . . . . . . . . . . 688
10.5 Eisenstat’s Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . 689
11 Miscellaneous Topics
690
11.1 Ordering for Coloring . . . . . . . . . . . . . . . . . . . . . . . . . . . 690
11.2 The Coloring Number . . . . . . . . . . . . . . . . . . . . . . . . . . . 694
11.3 Bounds for the Chromatic Number . . . . . . . . . . . . . . . . . . . . 696
11.4 Theoretical Results on Coloring Problems . . . . . . . . . . . . . . . . 696
12 Conclusion
698
12.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 698
12.2 Recent Development . . . . . . . . . . . . . . . . . . . . . . . . . . . . 699
12.3 Further Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 699

632

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

1. Introduction.
1.1. Overview. Many algorithms that solve nonlinear optimization problems anddiﬀerential equations require the computation of derivatives. For a vector function, thematrix of ﬁrst derivatives is the Jacobian; for a scalar function, the matrix of second
derivatives is the Hessian. These derivative matrices can be estimated through ﬁnite
diﬀerence approximations or computed exactly, within the limits of machine precision,
through automatic diﬀerentiation (AD). When the problems are large, many of theelements in the Jacobian or Hessian are zero, and these matrices are said to be sparse.
Sparsity and symmetry in the derivative matrices can be exploited to compute thenonzero entries eﬃciently. In the context of ﬁnite diﬀerence approximations, eﬃciency
corresponds to reducing the number of function evaluations required. In the context
of AD, eﬃciency often corresponds to reducing the number of AD passes required.
The problem of minimizing the number of function evaluations/AD passes needed
to compute a sparse derivative matrix can be formulated as a matrix partitioning
problem. Indeed, depending on the conditions and the methods employed, one canidentify ten variant matrix partitioning problems.
Each of these ten matrix partitioning problems can be modeled as a (specialized)
graph coloring problem. The graph coloring models provide deeper insight into thepartitioning problems and enable the design of eﬀective approximation or heuristic
algorithms. In this paper we revisit the research in this area of the last twenty
years with the aim of providing a unifying algorithmic framework. We integrate
existing coloring formulations with new ones, design more eﬃcient algorithms forproblems that have been studied earlier, and introduce graph models for new matrix
computation problems that arise in preconditioning.
1.2. A Jacobian Computation Problem. We begin the discussion of graph coloring models for the various matrix partitioning problems by considering a representative problem, a problem that occurs in a columnwise computation of a Jacobianmatrix. After providing a brief background of the problem within the context of aﬁnite diﬀerence approach, we will introduce a natural graph coloring formulation.
Consider a nonlinear function F : Rn → Rm , with n = m = 5. Let the nonzero
structure of the Jacobian matrix J of the function F be as shown in Figure 1.1(a).
In the ﬁgure, jk denotes the (k, ) entry of J when it is nonzero; the symbol “0”denotes a zero entry of J. (The illustration in Figure 1.1 is motivated by a ﬁgurefrom Hovland’s recent talk [67].)
Let ek ∈ Rn correspond to the kth coordinate vector (a vector with 1 in the kth
row and 0 in all other rows); then from the approximation 1 [F (x + ek ) − F (x)] ≈J(x)ek , by diﬀerencing the function along the direction ek , we can estimate the kth
column of J through the additional function evaluation F (x + ek ) (assuming F (x)has already been evaluated). Here is a small step size. A centered diﬀerence formula
would be more accurate, but for pedagogical purposes, the forward diﬀerence formulasuﬃces. Thus, if sparsity is not exploited, the estimation of a Jacobian matrix withn columns would require n additional function evaluations.
Now consider a subset of the columns of the Jacobian such that no two columns
have a nonzero in a common row; such a subset of columns is structurally orthogonal.
In a group of structurally orthogonal columns, the columns are pairwise orthogonal
to each other independent of the numerical values of the nonzeros. In the example,
columns 1 and 3 are structurally orthogonal; so are columns 1 and 4, and columns
3 and 5. Choose a column vector d with 1’s in components corresponding to theindices of columns in a structurally orthogonal group of columns, and zeros in all

633

GRAPH COLORING FOR COMPUTING DERIVATIVES
j 11 j 12 0

j 11 j 12

0

0

j 15

j 23 00000

j 32 j 33 j 34 0

j 41 00000

j 54 j 55

0

0

0

j 15

j 11 j 12 0

j 15

j 23 00

j 23 00000

j 32 j 33 j 34 0

j 41 00000

j 54 j 55

000

j 33 j 32 j 34 0
j 41 000

j 54 j 55

0

c1

c5

c2

0
0

0

0

j 15

j 11 j 12 j 15

j 23 000

j 32 j 33 j 34 0

j 41 0

0

0

0

0

0

0

j 54 j 55

0

j 41 00

j 54 0j 55

Fig. 1.1

b

c

c

3

c

4

c

5

c

r1

c1

r

2

c

3

c

4

c

5

c

r

c1

c

c2

5

r
r
r

c4

a

2

r

c3

j 23

j 34 j 32 j 33

c1

r

r

0

c4
j 11 j 12 0

r1

c3

d

2
3
4
5

2
3
4
5

e

Two partitions of a Jacobian into groups of structurally orthogonal columns (a)–(c). Eachpartition is also represented as a distance-1 coloring in the column intersection graph (d),and as a partial distance-2 coloring in the bipartite graph (e). In the upper row, four colors
are used: columns 1 and 3 are assigned the ﬁrst color, column 2 is assigned the second
color, column 4, the third color, and column 5, the fourth color. In the lower row, threecolors are used: columns 1 and 4 are assigned the ﬁrst color, column 2 is assigned thesecond color, and columns 3 and 5 are assigned the third color.

other components. The vector d is binary since each entry is either one or zero. By
diﬀerencing the function F along the vector d, one can simultaneously determine thenonzero elements in all of these columns through one additional function evaluation
at F (x + d).
Further, by partitioning the columns of the Jacobian into the fewest groups,
each consisting of structurally orthogonal columns, the number of (vector) function
evaluations needed to estimate the Jacobian matrix is minimized. This leads to thefollowing general problem: given the sparsity structure of a Jacobian matrix, partition
the columns into the fewest groups of structurally orthogonal columns.
Figure 1.1(b) shows two diﬀerent ways in which the Jacobian matrix can bepartitioned into groups of structurally orthogonal columns. In the ﬁgure, columnsof the same color belong to the same group. In the lower row of the ﬁgure, thematrix is partitioned into three such groups: {1, 4}, {2}, and {3, 5}. As we will see insection 3.5, the fewest groups of structurally orthogonal columns this matrix can bepartitioned into is three, the maximum number of nonzeros in a row.
Two “compressed” representations of the Jacobian matrix in which the nonzeros
in a group of structurally orthogonal columns are packed into the same column are
shown in Figure 1.1(c). Each column in the compressed representation of the Jacobian
corresponds to the nonzeros that can be simultaneously computed through one ﬁnite
diﬀerence operation.
Figure 1.2 shows a Jacobian matrix having 100 columns and its compressed representation with 28 columns obtained by grouping together nonzero elements in aset of structurally orthogonal columns into a single column. This ﬁgure is similar
to Figures 1.1(b) and (c); the nonzeros within each group of structurally orthogonal
columns are assigned a single color distinct from the colors assigned to other nonzeros, and then all nonzeros of the same color are grouped into one compressed column.
(Figure 1.2 is inspired by similar ﬁgures in Griewank [53] and Hovland [67].)

634

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

Fig. 1.2 A Jacobian matrix and its compressed representation.

Curtis, Powell, and Reid [35] were the ﬁrst to observe in 1974 that sparsity can be
employed in this way to reduce the number of function evaluations needed to estimate
the Jacobian. They described a simple greedy algorithm for partitioning the columnsof a Jacobian into a small number of groups of structurally orthogonal columns.
In 1983 Coleman and Moré [30] modeled this matrix partitioning problem as a
distance-1 graph coloring problem. The model uses the column intersection graphof a matrix where columns correspond to vertices and two vertices are joined by anedge whenever the corresponding columns have nonzeros in a common row (i.e., thecolumns are structurally nonorthogonal). A distance-1 coloring of a column intersection graph, i.e., an assignment of colors to the vertices such that adjacent (distance-1
neighboring) vertices receive diﬀerent colors, partitions the columns into groups ofstructurally orthogonal columns. Figure 1.1(d) shows the column intersection graphof the Jacobian in the example and illustrates two diﬀerent distance-1 colorings of thatgraph. The lower part of the ﬁgure shows that vertices c1 and c4 are assigned one color;vertex c2 is assigned a second color; and vertices c3 and c5 are assigned a third color.
The objective in all coloring problems considered in this paper, including the
distance-1 coloring problem, is to use the fewest colors possible. Since the distance1 graph coloring problem is known to be NP-hard, the work of Coleman and Moréshowed that it is highly unlikely that there is a polynomial time algorithm for partitioning the columns of a matrix into the fewest groups of structurally orthogonal
columns. Meanwhile, they developed several practically eﬀective heuristics for the
problem.
In this work, we use a diﬀerent graph coloring model for the same matrix partitioning problem. This coloring formulation uses a bipartite graph to represent a
Jacobian matrix. The vertex set V1 in the bipartite graph corresponds to the rows of
the matrix and the vertex set V2 corresponds to the columns. An edge joins a row
vertex rk to a column vertex c if the matrix element jk of the Jacobian is nonzero.

GRAPH COLORING FOR COMPUTING DERIVATIVES

635

In Figure 1.1, part (e) shows the bipartite graph of the Jacobian whose structure is
depicted in part (a).
Notice that two column vertices that are at a distance of two edges from eachother in the bipartite graph (hence they share a common row vertex as a neighbor)
are structurally nonorthogonal in the matrix. In fact, two columns are structurally
orthogonal if and only if they are at a distance greater than 2 from each other inthe corresponding bipartite graph. Thus, a distance-2 coloring of the set of column
vertices V2 is equivalent to a partitioning of the columns of the matrix into groupsof structurally orthogonal columns. A distance-2 coloring of the vertex set V2 is an
assignment of colors to these vertices such that every pair of column vertices at a
distance of exactly two edges from each other receives distinct colors. More precisely,
this coloring is a partial distance-2 coloring of the bipartite graph since the row vertex
set V1 is left uncolored. Figure 1.1(e) shows two diﬀerent distance-2 colorings of thecolumn vertices of the bipartite graph of the Jacobian in the example. These colorings
are the equivalents of the distance-1 colorings shown in part (d).
In the Jacobian estimation problem under discussion, the bipartite graph–based
partial distance-2 coloring formulation may be preferred over the column intersection
graph–based distance-1 coloring formulation for several reasons. One of the advantages of the former is that the distance-2 coloring problem serves as the most general
problem for the many variations of coloring problems addressed in this paper. The
identiﬁcation of such a general problem simpliﬁes the design of algorithms for the
specialized variants. Another advantage is related to the size of the graph used torepresent a matrix. The size of the column intersection graph of a matrix A is proportional to the number of nonzeros in AT A, whereas the size of the bipartite graphof A is proportional to the number of nonzeros in A. For many large-scale practical problems, the former could be considerably larger in size. A larger graph size inturn implies that more storage space is required and performance is degraded dueto the hierarchical structure of memory in modern computers [60]. In section 3.6,
we provide computational evidence to support this claim. A third advantage of the
bipartite graph–based formulation is related to ﬂexibility. Unlike the column intersection graph, the bipartite graph of a matrix is a faithful representation of the sparsity
structure and hence allows for a columnwise, rowwise, or combined column- and rowwise computation of a matrix. A detailed discussion of the comparison between the
two approaches is included in section 3, a section devoted to the formal treatment of
the problem under discussion. In section 9, we consider hypergraph coloring formulations for the same problem, where the equivalence between the two graph coloring
formulations can be viewed from yet another perspective.
1.3. Automatic Differentiation. For pedagogical reasons, the matrix partitioning problem introduced in section 1.2 has been discussed within the context of a ﬁnite
diﬀerence technique, a technique that computes an approximation to the derivative.
The same partitioning problem also arises when a Jacobian matrix is calculated using
the relatively recent technique of automatic diﬀerentiation (AD).
AD is a chain rule–based technique for evaluating the derivatives of functions
deﬁned by computer programs. Unlike ﬁnite diﬀerencing, the derivatives in AD arecomputed analytically, and hence without any truncation error.
AD has two basic modes of operation known as forward and reverse. These modes
correspond to a bottom-up and a top-down strategy of accumulating partial derivatives of elementary functions that deﬁne the computational scheme of the function to
be diﬀerentiated.

636

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

One pass of the forward mode of AD can be used to compute a column (or a groupof structurally orthogonal columns) of a derivative matrix and a pass of the reverse
mode can be used to compute a row (or a group of structurally orthogonal rows).
Hence, just as in the ﬁnite diﬀerence setting, the eﬃcient computation of a Jacobianmatrix using the forward mode of AD requires the columns to be partitioned into
as few groups of structurally orthogonal columns as possible. Similarly, an eﬃcient
application of the reverse mode requires a row-partitioning into the fewest possiblegroups of structurally orthogonal rows.
A number of AD software tools are available, including ADIFOR [16], Odyssee[110], and TAMC [50] for FORTRAN 77, and ADIC [17] and ADOL-C [55] for
C/C++.
Some AD tools, such as ADIFOR and ADIC, support a vector mode in which
several bidirectional derivatives are computed simultaneously. Using such tools, thecompressed derivative matrix (see Figure 1.2) can be computed in one pass. The computational eﬀort involved in these contexts is proportional to the number of columns
in the compressed matrix. Hence a partition into the fewest groups of structurally
orthogonal columns continues to capture the requirement of an eﬃcient computation.
1.4. Variations on Matrix Computation. As mentioned earlier, depending on
the type of derivative matrix being computed and the speciﬁcs of the method being
applied, there exist several variant matrix partitioning problems. Speciﬁcally, thenature of a particular problem in our context depends on the following:• whether the matrix to be computed is a Jacobian (nonsymmetric) or a Hessian (symmetric);• whether the evaluation scheme employed is direct or substitution-based (adirect method requires solving a diagonal system and a substitution method
relies on solving a triangular system of equations);• whether a unidirectional (1d) partition or a bidirectional (2d) partition isused (a unidirectional partition involves only columns or rows, whereas abidirectional one involves both columns and rows);• whether all of the nonzero entries of the matrix or only a subset need to be
determined; we refer to these as full and partial matrix computation.
Under certain assumptions, the aforementioned four mutually orthogonal factors
give rise to ten variants of matrix partitioning problems. Figure 1.3 depicts a classiﬁcation of the problems using these factors; the ﬁrst eight of the ten listed problems
are studied in this paper, the problems being numbered in the order in which they
appear in the body of this paper.
Each of these matrix partitioning problems can be modeled as a graph coloring
problem. In section 1.2 we have seen the graph coloring formulations of problem
P1, a partitioning problem arising in the computation of all nonzero entries of a
Jacobian matrix using a unidirectional partition via a direct method. In sections 1.4.1
through 1.4.4, we brieﬂy review the historical development of the coloring formulations
for problems P2 through P5 and introduce the new problems P6 through P8. The
discussion is organized along the four factors underlying the classiﬁcation in Figure 1.3.
Problems P1 through P8 are formally treated in detail in sections 3 through 8.
Note that the classiﬁcation in Figure 1.3 is oblivious as to whether the chosen
numerical technique is ﬁnite diﬀerence or AD. For example, as has been pointed
out earlier, problem P1 arises within the context of ﬁnite diﬀerence or AD. In our
discussion in what follows, we may refer to one numerical technique that simpliﬁes
the presentation, but the reader should bear in mind that the particular partitioning
problem also arises within the context of the other numerical technique.

637

GRAPH COLORING FOR COMPUTING DERIVATIVES

Full
2d
partition

1d partition

DirectPartialSubstitution

DirectSubstitution

2d
partition

1d partition

DirectSubstitution

DirectSubstitution

Jac.

Hes.

Hes.

Jac.

Jac.

Jac.

Hes.

Hes.

Jac.

Jac.

P1

P2

P4

P3

P5

P6

P7

P9

P8

P10

Fig. 1.3

Classiﬁcation of partitioning/coloring problems arising in sparse derivative matrix computation. Problems P1 through P8 are addressed in this paper.

1.4.1. Jacobian vs. Hessian. In 1979 Powell and Toint [107] extended the approach of Curtis, Powell, and Reid to compute sparse Hessians. Their approach
exploits the fact that only one of each pair of symmetric oﬀ-diagonal elements in the
Hessian matrix needs to be computed.
McCormick [97] introduced a distance-2 graph coloring model for the computation of Hessians in 1983. The model uses the adjacency graph representation of theunderlying symmetric matrix and requires that in every path u, v, w in the graph,vertices u, v, and w receive distinct colors. Independently, in 1984, Coleman andMoré [31] gave a more precise coloring model that exploits symmetry. They called
the model path coloring and required it to satisfy the two conditions: (1) every pairof adjacent vertices receives distinct colors (a distance-1 coloring), and (2) every pathon four vertices uses at least three colors. This variant of coloring was called starcoloring, e.g., by Fertin, Raspaud, and Reed [41], since in such a coloring every subgraph induced by vertices assigned any two colors is a collection of stars. We will use
the name star coloring in this paper since path coloring is a term used in the optical
networks literature to refer to a diﬀerent concept.
Unlike coloring problems corresponding to general matrices, coloring problems indiscretization of structured grids can be solved optimally in polynomial time. In 1984,
Goldfarb and Toint [52] studied such optimal estimation of Jacobians and Hessians
that arise as model problems in ﬁnite diﬀerence approximations of partial diﬀerential
equations.
1.4.2. Direct vs. Substitution-Based Evaluation. The evaluation scheme considered thus far is a direct one, in that each nonzero element of a derivative matrix is
obtained from some row of a matrix-vector product via a ﬁnite diﬀerence operationor an AD pass. Suppose A is an m × n derivative matrix and d1 , d2 , . . . , dp are binary
vectors (directions) corresponding to a column partition of A into p groups of structurally orthogonal columns. In the AD literature, the n × p matrix D consisting of
the p columns d1 , d2 , . . . , dp is known as a seed matrix. In a direct evaluation scheme,
the product AD can be used to deﬁne a diagonal system of equations where the unknowns are the nonzero entries of A and the right-hand side is the result obtained byp function evaluations or p passes of the forward mode of AD.
An alternative evaluation scheme would be substitution-based, in which the unknown matrix elements were determined by solving a triangular system of equations.

638

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

The latter is obtained by an appropriate (less restricted) choice of a seed matrix thatpartitions the columns of the derivative matrix into groups of not necessarily structurally orthogonal columns. There is a trade-oﬀ in the choice between a direct and asubstitution method. The former being more restrictive requires more function evaluations (or AD passes), while the latter is subject to numerical instability. Moreover,unlike a direct method where the unknowns are obtained without any further arithmetic, a substitution method incurs the additional computational cost involved withsolving a sparse triangular system of equations.
A substitution-based evaluation is often eﬀectively combined with the exploitationof symmetry and hence is used in computing the Hessian. A substitution method forcomputing the Hessian leads to an acyclic coloring problem in which the requirements
are that (1) the coloring corresponds to a distance-1 coloring, and (2) vertices in every
cycle of the graph are assigned at least three distinct colors. This variant of coloringis called acyclic since every subgraph induced by vertices assigned any two colors is acollection of trees—and hence is acyclic. Substitution methods for computing Hessians
were studied by Powell and Toint [107] and Coleman and Moré [31]. The latter authors
in particular found a less accurate model (in comparison with acyclic coloring) called
triangular coloring. The acyclic coloring formulation is due to Coleman and Cai [27],who called the coloring variant cyclic. The name acyclic coloring was introduced byGrünbaum [56], the ﬁrst to deﬁne and study the problem within an entirely diﬀerent
context.
A variant of a substitution method for computing a Jacobian matrix has been
considered by Hossain and Steihaug [63]. The method relies on ﬁrst ﬁnding a partitionof the columns into groups of structurally orthogonal columns and then merging apair of successive groups to get a column grouping that allows overlaps; in this way acolumn partition consisting of p groups is changed into a grouping having p−1 groups.
1.4.3. Unidirectional vs. Bidirectional Partition. The third source for problem
variation while computing the Jacobian lies in one’s choice to compute by partitioning
the columns (as we have discussed earlier), or rows, or both columns and rows. If
the Jacobian has a few dense columns, i.e., columns with a large number of nonzeros,
it may be more advantageous to consider a partitioning of the columns. Similarly, if
the matrix has a few dense rows, a partitioning of the rows may pay oﬀ. However,if the matrix contains a few dense columns and rows, it may be better to consider
partitioning subsets of both columns and rows. A partition that involves only columns
or rows is referred to as unidirectional, and one that involves both columns and rows
is called bidirectional. More speciﬁcally, given the structure of an m × n Jacobianmatrix A, a bidirectional partitioning problem is concerned with ﬁnding seed matrices
D1 ∈ {0, 1}n×p1 and D2 ∈ {0, 1}m×p2 such that the products AD1 and D2T A together
enable the determination of the entries of A and the value p = p1 + p2 is minimized.
Due to symmetry, there is no advantage in considering a bidirectional partition of the
Hessian; i.e., a symmetry-exploiting unidirectional partition suﬃces.
In the context of AD, bidirectional partitions arise when the Jacobian is computedby using the forward and reverse modes simultaneously.
Bidirectional partitioning of the Jacobian leads to specialized bicoloring problems in the bipartite graph, i.e., a coloring of subsets of both the row vertices andthe column vertices with disjoint sets of colors. When bidirectional partitioning isused within a direct evaluation scheme for Jacobians, the coloring problem is thatof star bicoloring; the corresponding model within a substitution-based scheme isthe acyclic bicoloring problem. Bidirectional partitioning problems and their graph

639

GRAPH COLORING FOR COMPUTING DERIVATIVES
Table 1.1

Graph coloring formulations for computing all nonzero entries of derivative matrices.
The Jacobian is represented by its bipartite graph, and the Hessian by its adjacency
graph. NA stands for not applicable. The labels in parentheses correspond to those usedin Figure 1.3.

Jacobian
HessianJacobian
Hessian

1d partition
Distance-2 coloring [P1]Star coloring [P2]NA
Acyclic coloring [P4]2d partition
Star bicoloring [P3]NA
Acyclic bicoloring [P5]NA

DirectDirectSubstitution
Substitution

coloring formulations were studied by Hossain and Steihaug [62] and Coleman andVerma [32].
1.4.4. Full vs. Partial Computation. The ﬁnal variation within the classiﬁcation
scheme of Figure 1.3 is whether all elements of the Jacobian and the Hessian are
required, or only a subset that is needed for preconditioning purposes. We refer tothese variations as full and partial matrix computation. The latter would be usefulin “matrix-free” methods for large-scale problems, where the Jacobian is too large
to be explicitly estimated, but a coarser representation of the Jacobian is used asa preconditioner. Partial matrix computation problems lead to restricted coloring
problems where only a speciﬁed subset of the vertices need to be colored; however,
one still needs to pay attention to the remaining vertices, since they could interferewith the estimation of the required matrix elements. Partial matrix computation
problems and their coloring formulations are treated in this paper for the ﬁrst time.
All of these variations lead to a rich collection of graph coloring problems. Table 1.1 shows the collection of ﬁve coloring problems that arise when we consider thecomputation of all nonzero entries of Jacobians and Hessians. Partial matrix computation problems lead to another set of ﬁve coloring problems, of which we haveformulated graph models only for direct methods, i.e., problems P6 through P8 in
Figure 1.3. The precise formulation of each of these coloring problems is described in
the section of the paper devoted to the corresponding problems.
1.4.5. Other Variations. One common feature in the partitioning/coloring problems cataloged in Figure 1.3 is that the nonzero entries of a derivative matrix are
obtained by solving either a diagonal or a triangular system of equations. A second
common feature is that the system of equations is produced from a partition of a setof columns or rows (or both) of a matrix; i.e., a set is divided into groups in which anelement of the set belongs to one and only one group. Newsam and Ramsdell [103]suggested a method for computing a Jacobian that relies on solving an overdeterminedrectangular system of equations, a system produced from a column grouping that is
not necessarily a partition. The approach of Newsam and Ramsdell has been applied
within the context of AD by Geitner, Utke, and Griewank [49]. Recently, Hossain
and Steihaug [64] suggested a method for computing a Jacobian that requires solving
a banded rectangular system of equations.
A third common feature in the problems of Figure 1.3 is that the granularity level
used in deﬁning the problems is an entire column or row of a matrix. In other words,
a column or a row is not divided any further. However, approaches that rely on more
“ﬁne-grained” partitioning have also been suggested. The element isolation methodof Newsam and Ramsdell [103] and the segmented column method of Hossain and
Steihaug [61, 65] are examples of such approaches.

640

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

The works mentioned in the previous two paragraphs are brieﬂy discussed in
section 10 of this article.
1.5. Objectives. Since the work of Curtis, Powell, and Reid [35], structural orthogonality in a sparse derivative matrix has been identiﬁed as a basic property to
be exploited for eﬃcient computation. Two columns of a Jacobian matrix are structurally orthogonal if and only if they are at a distance greater than two in the bipartite
graph of the Jacobian. Similarly, two columns of a Hessian matrix are structurally
orthogonal if and only if they are at a distance greater than two in the adjacency
graph of the Hessian. Thus distance-2 coloring of the bipartite graph of the Jacobianand the adjacency graph of the Hessian are the two archetypal problems in eﬃcientderivative matrix computation. This perspective emphasizes the central role playedby distance-2 coloring and thus provides a unifying algorithmic framework for the
many coloring problems.
Describing the uniﬁed algorithmic framework is the ﬁrst objective of this paper;
this framework simpliﬁes the development of algorithms and software for the problems
that have been considered in earlier work and clariﬁes the development of models andalgorithms for the new partial Jacobian and Hessian computation problems introduced
here. Many of the algorithms proposed earlier for the variant problems relied on
adding new edges to the graph until a distance-1 coloring of the transformed graphcould be used to solve the problem. We show that more eﬃcient algorithms are
possible for the variant problems by employing coloring algorithms directly to thenatural graph representations of the problem.
A second objective of this paper is to provide a self-contained introduction to theuse of graph coloring in derivative matrix computations, in order to make this work
accessible to a broad audience. For researchers in optimization, this article provides
the ﬁrst review of graph coloring methods for eﬃciently computing derivatives fornonlinear problems. The article has been organized by the manner in which thesecomputational methods would be used in optimization: Jacobian or Hessian matrix;
direct or substitution method; unidirectional or bidirectional estimation; full matrix
or partial matrix evaluation. Our hope is that this organization makes it easy forpractitioners to read the subset of topics of interest to them. An alternate way oforganizing this material would have been by the graph algorithms employed to solve
the problems, and indeed, that is how we ﬁrst attempted to present this material.
Work done in other research communities on graph coloring topics intersect withthe problems in derivative matrix computation. For instance, early results obtainedby the graph theory community on vertex orderings for reducing the number of colors
have been applied in the optimization context. The graph theory community has also
recently considered some of the variant coloring problems discussed here; e.g., starcoloring, which arises in Hessian estimation [5, 41, 102]. However, this community
seems to be unaware of work done nearly twenty years earlier on this problem in optimization! Conversely, earlier work on acyclic coloring was not noticed by researchersin optimization who discovered its application to evaluating derivative matrices with
substitution methods.
Graph and hypergraph coloring have been used in a wide collection of application areas in addition to optimization: register allocation in compilers [24], radio andwireless networks [44, 84], scientiﬁc computing [51, 73, 74, 111], data movement in distributed and parallel computing [58], facility location problems [117], cache-eﬃcientalgorithms [77], etc. Parallel computers make it feasible to solve large-scale problems
in many of these application areas, especially optimization, and hence there is cur-

GRAPH COLORING FOR COMPUTING DERIVATIVES

641

rently increased interest in eﬃcient algorithms and software for coloring graphs withmillions of vertices.
For the reasons spelled out in the previous two paragraphs, wider disseminationof this work among researchers in the optimization, scientiﬁc computing, computerscience, discrete mathematics, and various applications communities would be beneﬁcial.
A third objective of this work is to help dispel the feeling that most of the interesting work on this topic has been done already by pointing to the research frontier in thisarea: new coloring and ordering algorithms need to be developed and implemented formany problems considered in this article; new mathematical results that characterize
lower bounds for coloring special classes of graphs are needed. We believe that as thiswork becomes better known, new applications of coloring will be discovered, further
enriching the study of the mathematics and algorithms in this area.
1.6. New Contributions. This article gives a coherent review of earlier work, but
it also has several new results. The following summary highlights our contributions.
1. We propose distance-2 coloring as a natural, ﬂexible, space-eﬃcient model for
problem P1 (section 3).
2. We expose the inter-relationships among the coloring formulations of the various matrix partitioning problems and identify distance-2 coloring as a unifying, archetypal model (section 7, building on the discussions in sections 3–6).
3. Using the insight gained from the uniﬁed graph-theoretic treatment, we develop several greedy heuristic algorithms, including new algorithms for starcoloring (section 4) and star bicoloring (section 5). We also sketch a newacyclic coloring algorithm that exploits the structure of two-colored induced
subgraphs (section 6). In a separate paper (see [48] for a draft), we will describe this algorithm in more detail, as well as a new star coloring algorithm
based on the same paradigm.
4. We report experimental results from our implementations that demonstrate• the advantages oﬀered by the bipartite graph–based distance-2 coloring
formulation of problem P1 in comparison with the column intersection
graph–based distance-1 coloring formulation (section 3), and• the time/quality trade-oﬀ between two star coloring algorithms; one of
the algorithms is new while the other is a translation from a previously
known matrix-based algorithm (section 4).
5. In the context of bicoloring models, we make the connection to the problem
of ﬁnding a vertex cover in a graph explicit (section 5).
6. We develop the ﬁrst graph coloring formulations for partial matrix computation problems, problems where only a subset of the nonzero entries of amatrix is required to be computed for preconditioning purposes (section 8).
7. We give the ﬁrst hypergraph coloring formulations for problem P1. These
formulations provide an interesting alternative perspective and are likely to be
useful in the context of partial Jacobian and Hessian computation (section 9).
8. We formulate the element isolation technique of Newsam and Ramsdell [103]for computing a Jacobian as an edge coloring problem in the associated bipartite graph (section 10).
1.7. Scope and Further Reading.
Automatic Differentiation. A discussion of the technical details of AD is beyond thescope of this paper. The book by Nocedal and Wright [104] contains an introductory
discussion of AD and how graph coloring can be used to reduce the number of passes

642

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

Table 1.2 A list of concepts known by diﬀerent names. The names acyclic and star coloring were
earlier used by Grünbaum [56] and Fertin, Raspaud, and Reed [41], respectively.

HereStar coloring

Acyclic coloring
Star bicoloring
Acyclic bicoloring
Structurally orthogonal partition
Symmetrically orthogonal partition
Bidirectional partition

Elsewhere
Symmetric coloring [31]Path coloring [27]Distance- 32 coloring [45, 47]
Cyclic coloring [27]Path bicoloring [32]Cyclic bicoloring [32]Consistent partition [30]Symmetrically consistent partition [31]Bipartition [32]Row-column partition [62]needed to compute the Jacobian and the Hessian. The recent book by Griewank [53]has an excellent discussion of the principles and techniques of AD. The books [33,54] contain several papers that discuss diﬀerent aspects of AD, including its theory,
implementations, and applications. Examples of works that discuss AD within thecontext of parallel computation include [68, 69, 106].
Combinatorial problems other than coloring also occur in optimizing the computation of derivatives using AD. A graph model for eliminating vertices and edges in acomputational graph can be used in this context to reduce the storage and the eﬀort
needed to compute the Jacobian or the Hessian. Depending on the type of elimination
technique employed, there exist a number of combinatorial problems that could be
regarded as searching for shortest paths in graphs. These problems are similar tograph models of sparse Gaussian elimination. Naumann [101] has recently studied
graph elimination in AD. Examples of studies that investigate such methods can also
be found in [33, 53, 54, 100]. In a recent talk, Hovland [67] surveyed the various
combinatorial problems that occur in AD. We will not discuss these problems in thispaper.
Graph Coloring. The graph coloring literature per se is vast. This work focuseson coloring problems that arise in computing derivative matrices. Furthermore, ouremphasis is on problem formulations, graph models, and fast heuristic or approximation algorithms. Iterative local improvement heuristics and exact algorithms for thecoloring problems are not discussed. In section 11 we include a discussion of eﬀective
vertex orderings for greedy distance-1 coloring, a discussion of the coloring numberof a graph, and a brief review of theoretical results on some of the coloring problems
addressed in this paper. The monograph by Jensen and Toft [71] provides a detailed
account of over 200 open coloring problems and an extensive list of references. The
survey article by Toft [115] has a good discussion of mathematical results for the
distance-1 coloring problem. The recent book edited by Kubale [86] includes survey
articles on several diﬀerent coloring problems. The WWW pages of Culberson [34]and Trick [116] are two useful web resources for the distance-1 coloring problem.
1.8. Terminology. For some of the concepts described in this paper, we haveintroduced names that diﬀer from those used earlier in the numerical optimization
literature. We have done so when we felt that the new names are more descriptive or
consistent with the terminology in other ﬁelds. For easy reference, Table 1.2 lists theimportant concepts in this paper that have been named diﬀerently elsewhere.
2. Deﬁnitions of Graph-Theoretic Concepts. In the current section we deﬁne
most of the graph-theoretic concepts used throughout this paper. Other concepts will
be deﬁned later as required.

GRAPH COLORING FOR COMPUTING DERIVATIVES

643

2.1. Preliminary Concepts and Notations. A graph G is an ordered pair (V, E),
where V is a ﬁnite and nonempty set of vertices and E is a set of unordered pairs of
distinct vertices called edges. If (u, v) ∈ E, vertices u and v are said to be adjacent;otherwise they are called nonadjacent. A path of length (edges) in a graph is asequence v1 , v2 , . . . , v+1 of distinct vertices such that vi is adjacent to vi+1 for 1 ≤i ≤ . Two distinct vertices are said to be distance-k neighbors if a shortest pathconnecting them has length at most k. Note that the set of distance-k neighbors ofvertex u, denoted by Nk (u), does not include u itself. Note also that, by deﬁnition,
two distance-k neighboring vertices are also distance-k  neighbors for k  > k. The
usual deﬁnition of the degree of a vertex is the number of edges incident on it. We
extend this notion and deﬁne the degree-k of a vertex u, denoted by dk (u), to be thenumber of distinct paths of length at most k edges starting at u. The terms degree anddegree-1 are thus synonymous, and they are used interchangeably in this paper. The
average degree-k in a graph is denoted by δ k . For brevity, we often write d(u) insteadof d1 (u), and δ instead of δ 1 . The maximum degree-1 in a graph is (consistently)
denoted by ∆, and the minimum degree-1 is denoted by δ.
For a simple example consider the graph shown on the left-hand side of Figure 2.1.
In this example, the distance-1 and distance-2 neighbors of vertex a are N1 (a) = {b, c}and N2 (a) = {b, c, d, e}; the degree of vertex c is d1 (c) = 4, which is also the maximum
degree, ∆, in the graph. The degree-2 of the vertex c is d2 (c) = 8.
A graph is called bipartite if its vertex set can be partitioned into two disjoint
sets V1 and V2 such that every edge of the graph connects a vertex in V1 with avertex in V2 . In a bipartite graph, we denote the maximum degree in the vertex set
V1 by ∆(V1 ), and the maximum degree in the vertex set V2 by ∆(V2 ). Similarly, theaverage degree-k in the sets V1 and V2 are denoted by δ k (V1 ) and δ k (V2 ), respectively.
The graph in Figure 1.1(e) is an example of a bipartite graph; V1 consists of therow-vertices and V2 consists of the column-vertices. In this graph, ∆(V1 ) = 3 and∆(V2 ) = 2.
A graph is connected if any two of its vertices are linked by a path in the graph.
A tree is a connected graph without any cycle. A forest is a graph each of whose
components is a tree. A bipartite graph Gb = (V1 , V2 , E) in which each vertex in V1
is connected to every vertex in V2 is complete. A complete bipartite graph in which
one of the vertex sets consists of a single vertex is a star.
Let G = (V, E) be a graph. A set of vertices C ⊆ V is said to cover a set of edges
F ⊆ E if for every edge e ∈ F , at least one of the endpoints of e is in C. If the vertex
set C covers the set of edges E, it is called a vertex cover. A set of vertices I ⊆ V iscalled an independent set if no two vertices in I are adjacent to each other. A set ofvertices Q ⊆ V is called a clique if the vertices in Q are pairwise mutually adjacent
to each other. For any set of vertices U ⊆ V , the graph induced by U is the subgraphof G whose vertex set is U and whose edges are precisely the edges of G with both
endpoints in U . The graph induced by U is denoted by G[U ].
In the graph on the left in Figure 2.1, the set {a, c, f } is a vertex cover, the set{a, d, e} is an independent set, and the set {a, b, c} is a clique.
2.2. Distance-k Graph Coloring. A distance-k (vertex) coloring of a graph G =(V, E) is a mapping φ : V → {1, 2, . . . , p} such that φ(u) = φ(v) whenever vertices
u and v are distance-k neighbors. The color assignment to the vertices of the graphon the left in Figure 2.1 shows an example of a distance-2 coloring. The least possible number of colors required for a distance-k coloring of a graph G is called its
k-chromatic number and is denoted by χk (G). A distance-k coloring of G = (V, E)
is called partial if only a subset of the vertices is colored; in particular, a partial

644

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

a 1

c
3

d 4

Fig. 2.1

2 b

a 1

5 e

d 4

c
3

2 b

5 e

1

1

f

f

An example of a graph G (left) and its square graph G2 (right). The color assignments(numbers within the circles) show a distance-2 coloring of G and the equivalent distance-1
coloring of G2 .

distance-k coloring of G = (V, E) on W , W ⊂ V , is a mapping φ : W → {1, 2, . . . , p}such that φ(u) = φ(v) whenever vertices u and v in the set W are distance-k neighbors. Figure 1.1(e) shows two examples of a partial distance-2 coloring of a bipartite
graph on one of its vertex sets.
2.3. The Power of a Graph. The kth power of a graph G is the graph Gk whose
vertex set is the same as that of G and whose edge set consists of pairs of vertices(u, v) whenever vertices u and v are distance-k neighbors in G. For example, thegraph on the right in Figure 2.1 is the square graph G2 of the graph G shown on theleft. The following equivalence is easy to see.
Lemma 2.1. A mapping φ is a distance-k coloring of G if and only if it is adistance-1 coloring of Gk .
Figure 2.1 illustrates the equivalence between a distance-2 coloring of a graph anda distance-1 coloring of its square.
2.4. Representing Matrices Using Graphs. Here we will deﬁne diﬀerent graphrepresentations of the sparsity structure of matrices. First, a remark on notation isin order. In the rest of this paper, for a given matrix A, the ith row is denoted by ri ,the jth column is denoted by aj , and the (i, j) entry is denoted by aij . Moreover, inorder to simplify notation, we use ri (aj ) to refer both to the ith row (jth column)of matrix A and to the corresponding vertex in an appropriate graph representationof A.
Bipartite Graph. Let A be an m × n matrix with rows r1 , r2 , . . . , rm and columns
a1 , a2 , . . . , an . The bipartite graph Gb (A) of A is deﬁned as Gb (A) = (V1 , V2 , E),
where V1 = {r1 , r2 , . . . , rm }, V2 = {a1 , a2 , . . . , an }, and (ri , aj ) ∈ E whenever aij isnonzero for 1 ≤ i ≤ m, 1 ≤ j ≤ n. For example, Figure 1.1(e) shows the bipartite
graph representation of the matrix whose nonzero structure is shown in part (a).
The size of the graph Gb (A) is proportional to the size of the matrix A: thenumber of vertices |V1 | + |V2 | = m + n, and the number of edges |E| = nnz(A), wherennz(A) is the number of nonzeros in A.
Adjacency Graph. Let A be an n × n symmetric matrix having nonzero diagonal
elements with its columns denoted a1 , a2 , . . . , an . The adjacency graph of A is G(A) =(V, E), where V = {a1 , a2 , . . . , an }, and (ai , aj ) ∈ E whenever aij is nonzero for1 ≤ i ≤ n, 1 ≤ j ≤ n, i = j.
The graph G(A) exploits the symmetry available in A since both nonzero entries
aij and aji are represented by the single edge (ai , aj ). In the adjacency graph G(A),GRAPH COLORING FOR COMPUTING DERIVATIVES

645

the edges corresponding to the nonzero diagonal elements of A are not explicitly
represented. Thus, the number of edges in G(A) is 21 (nnz(A) − n), where nnz(A) isthe number of nonzeros in A. Since there is a vertex corresponding to each column,
the number of vertices in G(A) is n.
Column Intersection Graph. The graph formulations for the derivative matrix computation problems considered in this paper are based solely on the two graph representations discussed above. However, the column intersection graph has also been used
to represent a nonsymmetric matrix.
Let a1 , a2 , . . . , an now correspond to the structures of the columns of an m × n
matrix A. In particular, let each aj = {i ∈ {1, . . . , m} : aij = 0} (i.e., aj is a set of row
indices in which the jth column of A has a nonzero entry). The column intersection
graph of A is Gc = (V, E), where V = {a1 , . . . , an } and (ai , aj ) ∈ E ⇐⇒ ai aj = ∅.
In other words, an edge (ai , aj ) exists whenever both the ith and the jth column of A
have nonzero entries in at least one common row index. For example, Figure 1.1(d)shows the column intersection graph of the matrix whose nonzero structure is shownin part (a).
Note that the column intersection graph of a matrix A is isomorphic to the adjacency graph of AT A.
3. Unidirectional, Direct Computation of the Jacobian. In this section we revisit the partitioning problem introduced in section 1.2 with more technical details.
Given a continuously diﬀerentiable function F : Rn → Rm , the Jacobian of∂fi
F at the point x is the m × n matrix whose (i, j) entry F  (x)ij = ∂x(x), wherejf1 (x), f2 (x), . . . , fm (x) are the components of F (x). Let A denote the Jacobian matrix
F  (x). An estimate for the jth column of A can be obtained from the ﬁnite diﬀerence
approximation(3.1)

Aej = aj =∂1
F (x) ≈ [F (x + ej ) − F (x)],∂xj

1 ≤ j ≤ n,

where ej is the jth coordinate vector and is a positive step length. Formula (3.1)
shows that if each column of A is computed independently, n function evaluations,in addition to the evaluation of F (x), will be required. However, by exploiting thesparsity structure of A, the required number of function evaluations can be reduced
signiﬁcantly. Often the sparsity structure of A is either available or can be obtained
relatively easily [53, 61]. The goal here is to exploit the known sparsity structure of
A to estimate its nonzero entries using as few function evaluations as possible. We
assume that evaluating the function F at a given point is more eﬃcient than evaluating
the components fi , for 1 ≤ i ≤ m, separately. This assumption is reasonable sincein many applications the components of F have common subexpressions, and hence
the computation can be made eﬃcient by avoiding repeated computation of thesesubexpressions.
A sparsity-exploiting estimation of a matrix using ﬁnite diﬀerences involves ﬁnding directions (vectors) d1 , d2 , . . . , dp such that the products Ad1 , Ad2 , . . . , Adp enable
the determination of all the nonzero entries of A. Here each vector d is binary (i.e.,each entry is either 1 or 0) and is obtained via some linear combination of the coordinate vectors. Using (3.1), the products Ad1 , Ad2 , . . . , Adp can be used to deﬁne asystem of linear equations in which the unknowns are the nonzero elements of A. If
the choice of the binary vectors di is such that the resulting system of equations canbe ordered to be diagonal, then A can be determined directly. If the vectors di arechosen such that the system of equations can be ordered to be triangular, then theunknowns can be determined via substitution.

646

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

In the current and the next two sections, we consider direct methods; substitution
methods will be discussed in section 6.
3.1. The Matrix Partitioning Problem. In a direct method for determining a
Jacobian matrix A, for each nonzero element aij there should exist some vector dk inthe set {d1 , d2 , . . . , dp } such that aij = (Adk )i . Here (Adk )i denotes the ith entry of
the vector Adk . Thus each nonzero matrix element aij can be read oﬀ from an entry
of some vector Adk .
In an eﬃcient estimation of a Jacobian, the number of vectors p, which is proportional to the number of function evaluations required, needs to be minimized. Hencethe following general problem can be stated.
Problem 3.1. Given the sparsity structure of an m × n matrix A, ﬁnd the fewest
binary vectors d1 , d2 , . . . , dp such that the products Ad1 , Ad2 , . . . , Adp enable a direct
determination of A.
Curtis, Powell, and Reid [35] were the ﬁrst to address Problem 3.1. They observedthat the vectors d1 , d2 , . . . , dp can be obtained by partitioning the columns of A into
groups of structurally orthogonal columns.
Definition 3.2. A partition of the columns of a matrix A is structurally orthogonal if for every nonzero element aij the group containing column aj has noother column with a nonzero in row ri .
Let {C1 , C2 , . . . , Cp } be a structurally orthogonal column partition. With eachgroup Ck , associate a binary vector dk having components δj = 1 if aj belongs to Ck ,and δj = 0 otherwise. Then
Adk =

n

j=1

δj aj =

aj .

aj ∈Ck

If aij = 0 and column aj ∈ Ck , then aij = (Adk )i . Thus, all the nonzero entries of A
can be determined from the p matrix-vector products Ad1 , . . . , Adp .
Problem 3.1 can thus be cast as a matrix partitioning problem in the followingway.
Problem 3.3. Given the sparsity structure of an m × n matrix A, ﬁnd a structurally orthogonal partition of its columns that has the fewest groups.
As discussed in the introduction, Problem 3.3 also arises while computing thederivative matrix A using AD. In such a context, the number of groups in a structurally
orthogonal partition corresponds to the number of passes of the forward mode of AD.
3.2. A Graph Coloring Formulation. Recall that a distance-k coloring of agraph G = (V, E) is a mapping φ : V → {1, 2, . . . , p} such that φ(u) = φ(v) wheneveru and v are distance-k neighbors. The distance-k graph coloring problem asks for a
distance-k coloring with the fewest colors.
A distance-k coloring of G = (V, E) that uses p colors partitions the vertex set V
into p color classes U1 , U2 , . . . , Up , where Ui = {u ∈ V : φ(u) = i}. Each color classis a distance-k independent set; i.e., no pair of distinct vertices in the class consistsof distance-k neighbors. This prompts the question, Is there a natural relationshipbetween Problem 3.3 and the distance-k graph coloring problem? The following simple
observation gives a graph-theoretic characterization of structural orthogonality in anonsymmetric matrix, which leads to an answer to this question.
Lemma 3.4. Let A be an m × n matrix and Gb (A) = (V1 , V2 , E) be its bipartite
graph. Two columns (or rows) in A are structurally orthogonal if and only if thecorresponding vertices in Gb (A) are at a distance greater than two from each other.

GRAPH COLORING FOR COMPUTING DERIVATIVES

647

Proof. We prove the statement for columns; a similar argument can be used toprove the case for rows. First notice that in the bipartite graph Gb (A), any path
between a pair of distinct column vertices is of even length.
Assume that vertices ai and aj in V2 are at a distance greater than two fromeach other. This implies that there exists no path ai , 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 thatboth aki and akj are nonzero. Hence, by deﬁnition, columns ai and aj are structurally
orthogonal.
Conversely, assume that columns ai and aj are structurally orthogonal. Then, by
deﬁnition, there is no k ∈ [1, m] such that aki = 0 and akj = 0. This implies that
there is no path ai , rk , aj in Gb (A) for any 1 ≤ k ≤ m. Hence, vertices ai and aj areat a distance d > 2 from each other. In particular, vertices ai and aj are a distanced ≥ 4 apart.
By Lemma 3.4, ﬁnding a structurally orthogonal partition of the columns of amatrix A is equivalent to ﬁnding a partial distance-2 coloring of Gb (A) = (V1 , V2 , E)
on V2 . (The coloring is partial since V1 is not colored.) Theorem 3.5 states thisequivalence.
Theorem 3.5. Let A be an m × n matrix and Gb (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 structurally orthogonal partition of the columns of A.
In view of Theorem 3.5, Problem 3.3 is equivalent to the following graph coloring
problem.
Problem 3.6. Given the bipartite graph Gb (A) = (V1 , V2 , E) representing thesparsity structure of an m × n matrix A, ﬁnd a partial distance-2 coloring of Gb (A)on V2 that uses the fewest colors.
Note that for matrices with a few dense rows, a row partition may yield fewergroups than a column partition. Consequently, the matrix problem one needs to solve
is Problem 3.3 applied to AT . In such cases, the bipartite graph formulation becomes
handy—the equivalent problem is to ﬁnd a partial distance-2 coloring in the samegraph, but now on the vertex set V1 .
3.3. An Alternative Coloring Formulation. Coleman and Moré [30] were theﬁrst to formulate Problem 3.3 as a graph coloring problem. They showed Problem 3.3
to be equivalent to the distance-1 graph coloring problem on the column intersection
graph of the underlying matrix. Here we will show the equivalence between the
bipartite graph–based distance-2 coloring formulation and the column intersection
graph–based distance-1 coloring formulation of Problem 3.3.
The power of a graph gives an alternative view to the distance-k graph coloring
problem. Recall that the kth power of a graph G = (V, E) is the graph Gk = (V, F )where (u, v) ∈ F if and only if u and v are distance-k neighbors in G. As stated in
Lemma 2.1, a distance-k coloring of G is equivalent to a distance-1 coloring of Gk .
A particular implication of Lemma 2.1 is that a distance-2 coloring of a graph isequivalent to a distance-1 coloring of the square of the graph. This establishes theequivalence between distance-2 coloring the column vertices of the bipartite graphand distance-1 coloring the vertices of the column intersection graph of a matrix.
Speciﬁcally, the column intersection graph Gc (A) of a matrix A is isomorphic to the2adjacency graph of AT A. We note that Gc (A) is in fact the subgraph of Gb (A)induced by the vertices in V2 .
Lemma 3.7. Let Gb (A) = (V1 , V2 , E) and Gc (A) = (V2 , E  ) be the bipartite andcolumn intersection graphs of matrix A. Then Gc = Gb 2 [V2 ].

648

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN
a1
a1
r1

1

a2

1

1
a7

a3
a11 0
0
0
0 a16 a17 a18
0 a22 0
0 a25 0 a27 a28
0
0 a33 0 a35 a36 0 a38
00
0 a44 a45 a46 a47 0

r2

r3

a4
a5

1

2
1

a4

1

a6
1

1

2
2
2

2

a8
2

a2
1

a6
a5
r4

a7
a8

1

1

1

a3

Fig. 3.1 A matrix and its bipartite and column intersection graph representations.

3.4. Comparison. Figure 3.1 depicts a matrix A, the corresponding bipartite
graph Gb (A) = (V1 , V2 , E), and the column intersection graph Gc (A) = (V2 , E  ).
Note that we have augmented the graph Gc (A) with edge weights—w(ai , aj ) is thesize of the intersection of the sets represented by vertices ai and aj . In terms of thematrix A, the weight w(ai , aj ) is the total number of rows where both the ith andthe jth columns of A have nonzero entries. In the bipartite graph Gb (A), the weightw(ai , aj ) corresponds to the number of common neighbors of vertices ai and aj .
In sections 3.4.1 through 3.4.5, we compare and contrast the two graph formulations in terms of ﬂexibility, suitability for uniﬁcation, graph size, ease of graphconstruction, and use of existing software.
3.4.1. Flexibility. Notice that the column intersection graph is not an equivalent
representation of the sparsity structure of the underlying matrix. In particular, given
an edge between two column vertices, one cannot determine the row index (indices)
at which the columns share nonzero entries. In contrast, the bipartite graph is anequivalent representation of the structure of the matrix. This provides ﬂexibility. For
instance, the bipartite graph can be used in a column-only, row-only, or combined
row and column partition of the matrix. The column intersection graph, on the other
hand, is applicable only to a column partition. In general, the advantage of the
bipartite graph representation (in the context of computing derivative matrices) isthat the representation is decoupled from the eventual technique to be employed andthe matrix entries to be determined.
3.4.2. Uniﬁcation. As has been discussed in the introduction, depending on thetype of matrix to be evaluated and on the method employed, there exist several
variants of coloring problems in the context of derivative matrix computation. One of
the key advantages of the bipartite graph–based formulation for Jacobian estimationin this regard is that it leads to the identiﬁcation of an archetypal problem: thedistance-2 coloring problem. The identiﬁcation of an archetype is particularly useful
in developing new algorithms and software for the specialized variants.
3.4.3. Graph Size. Although Lemma 3.7 correlates the bipartite graph of a matrix with its column intersection graph, one cannot immediately deduce that one
graph is larger in size than the other. The size of the particular graph depends on thesparsity structure of the matrix. Here we provide a rough analysis to show that forsparse matrices of practical interest, the column intersection graph is likely to have
many more edges than the bipartite graph representation.

GRAPH COLORING FOR COMPUTING DERIVATIVES

649

Let A be a matrix, Gb (A) = (V1 , V2 , E) be the bipartite graph of A, and Gc (A) =(V2 , E  ) be the (weighted) column intersection graph of A. Further, let w denote theaverage edge weight in Gc . Thenw(e) =

e∈E 

1 
·
2



(d1 (v) − 1).

u∈V2 v∈N1 (u)Expressing the sum on the left-hand side with average weight and bounding d1 (v) by∆(V1 ), we get|E  | · w ≤1d1 (u).
(∆(V1 ) − 1)2u∈V2

Expressing the sum on the right-hand side with average degree gives|E  | · w = O(|V2 | · δ 1 (V2 ) · ∆(V1 )).
Since |V2 | · δ 1 (V2 ) = |E|, we have
(3.2)

|E  | = O(|E| · ∆(V1 )/w).

Therefore, since the inequality ∆(V1 ) > w is likely to hold in practice, and the constant in the Big Oh is close to 1, the column intersection graph of the matrix could
have many more edges than the bipartite graph. The larger graph size has in turntwo implications. First, more storage space would be required. Second, the averageruntime per edge spent in constructing the graph is likely to be larger since more data
has to be stored in memory, and memory access time is nonuniform on a hierarchicalmemory processor.
3.4.4. Ease of Construction. The sparsity structure of the matrix A is identicalto the adjacency lists of the vertices of the bipartite graph Gb (A), and hence theconstruction of the latter is trivial. In principle, the data structure used to represent
A could be used for implementing algorithms that use Gb (A). In contrast, Gc (A)has to be computed. The time required for the computation of Gc (A) is proportional
to the number of edges in Gc (A). Using (3.2) and noting that w ≥ 1, the following
lemma gives the complexity of constructing Gc (A).
Lemma 3.8. Given a graph Gb (A) = (V1 , V2 , E), the time required for constructing Gc (A) is O(|E| · ∆(V1 )).
It should, however, be noted that once the graph Gc (A) is computed, a subsequentdistance-1 coloring of Gc (A) can be done faster than a distance-2 coloring of Gb (A).
In section 3.5, we show that the overall time complexity of constructing and thendistance-1 coloring Gc (A) is of the same order as that of distance-2 coloring Gb (A).
3.4.5. Use of Existing Software. Serial program packages that implement various practically eﬀective distance-1 coloring heuristics exist [28, 29]. For matrix partitioning problems where a column intersection graph–based formulation can be applied,these packages can be readily used. On the other hand, since distance-2 coloring is anarchetypal model in our context, eﬃcient programs, including parallel ones, for the
distance-2 coloring problem need to be developed.
3.5. Algorithms. For every ﬁxed integer k ≥ 1, the distance-k graph coloring
problem is NP-hard [89]. Thus, in practice, one is bound to rely on approximation
algorithms or heuristics. An algorithm A is said to be a γ-approximation algorithm for650

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

Algorithm 3.1. A greedy distance-2 coloring algorithm.
procedure D2ColoringAlg(G = (V, E))
Let v1 , v2 , . . . , v|V | be a given ordering of V
Initialize forbiddenColors with some value a ∈ V
for i ← 1 to |V | dofor each colored vertex w ∈ N1 (vi ) do
forbiddenColors[color[w]] ← vi
for each colored vertex x ∈ N1 (w) do
forbiddenColors[color[x]] ← vi
end forend forcolor[vi ] ← min{c > 0 : forbiddenColors[c] = vi }end forend procedurea minimization problem if its runtime is polynomial in the input size and if for every
problem instance I with the value OP T (I) of an optimal solution, the value A(I) of
A(I)the solution produced by A is such that OP
T (I) ≤ γ. For a minimization problem,the approximation ratio γ satisﬁes γ ≥ 1, and the goal in designing an approximation
algorithm 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.
The current best known approximation ratio for the distance-1 coloring problem
is O(|V |( logloglog|V|V| | )3 ) [59]. Also, the problem is known to be not approximable within

O(|V |1/7− ) for any > 0, unless P = NP [12]. Despite these rather pessimistic
results, there exist several practically eﬀective distance-1 coloring heuristics [30, 94].
A detailed discussion of some of these heuristics is included in section 11 of thispaper. In the current section we show how a distance-1 coloring heuristic can be easily
adapted to the distance-2 coloring case by looking at the extended neighborhood of avertex. The algorithms we present in this paper are greedy in nature; i.e., the verticesof a graph are processed in some order and at each step a decision that looks best at
the moment (and that will not be changed later) is made.
In section 3.5.1 we present a generic greedy distance-2 coloring algorithm and givea detailed analysis of its performance in terms of both computation time and numberof colors used. In later sections, adaptations of this algorithm tailored to the various
coloring problems that concern us will be presented.
3.5.1. A Distance-2 Coloring Algorithm. A simple approach for an approximate
distance-2 coloring of a graph G = (V, E) is to visit the vertices in some order, ateach step assigning a vertex the smallest color that is not used by any of its distance-2
neighbors.
The degree-2 of a vertex v in G is bounded by ∆2 ; i.e., d2 (v) = w∈N1 (v) d1 (w) ≤∆ · d1 (v) ≤ ∆2 . The vertices in G can always be distance-2 colored trivially using |V |diﬀerent colors. Thus, it is always possible to assign a vertex a color chosen from theset {1, 2, . . . , min{∆2 + 1, |V |}}. D2ColoringAlg, outlined in Algorithm 3.1, uses
this fact as it colors the vertices of the graph in a given order. In the algorithm, color
is a vertex-indexed array that stores the color of each vertex and forbiddenColors
is a color-indexed array of size Cmax = min{∆2 + 1, |V |} used to mark the colors

GRAPH COLORING FOR COMPUTING DERIVATIVES

651

that cannot be assigned to a particular vertex. Speciﬁcally, forbiddenColors[c] = vi
indicates that color c cannot be assigned to vertex vi .
We prove the following statement regarding D2ColoringAlg.
Lemma 3.9. D2ColoringAlg ﬁnds a distance-2 coloring of a graph G = (V, E)
in time O(|V |δ 2 ).
Proof. We ﬁrst show correctness. In step i of the algorithm, the colors used byeach of the colored distance-2 neighbors of vertex vi are marked (using vi ) in the arrayforbiddenColors. Thus, at the end of the middle for-loop, the colors that are allowed
for vertex vi are the indices in forbiddenColors, where the marker used is diﬀerent
from vi . The minimum value in this set is thus the smallest allowable color for vertex
vi . Notice that the array forbiddenColors does not need to be initialized at every stepas the marker vi is used only in step i.
Turning to complexity, marking the forbidden colors at step i of the algorithm
takes O(d2 (vi )) time. Finding the smallest allowable color to assign vi can be done
within the same order of time by scanning the array forbiddenColors sequentially
until the ﬁrst index
 c where a value other than vi is found. The total time is thus
proportional to vi ∈V d2 (vi ) = O(|V |δ 2 ).
We now analyze the quality of the output of D2ColoringAlg. The distance-1
neighbors of a vertex in a graph G form a clique in the square graph G2 . This fact
immediately provides a lower bound on χ2 (G), the distance-2 chromatic number of G.
Lemma 3.10. For every graph G, χ2 (G) ≥ ∆ + 1.
In general, the graph G2 may contain a clique whose size is larger than ∆ + 1.
Thus, the maximum clique size in G2 is a tighter lower bound on χ2 (G) than thebound ∆ + 1.
Let the number of colors used by D2ColoringAlg on a graph G = (V, E) beχ2 (G, greedy). As discussed earlier, χ2 (G, greedy) ≤ min{∆2 + 1, |V |}. Combining
this with the lower bound on χ2 given in Lemma 3.10, we get the following theorem
and its corollary. The latter is due to McCormick [97].
2
Theorem 3.11. ∆ + 1 ≤ χ2 (G) ≤ χ2 (G, greedy) ≤ min{∆ + 1, |V |}. 
Corollary 3.12. D2ColoringAlg is an O(min{∆, |V |}) = O( |V |) approximation algorithm.
1
Proof. From Theorem 3.11, the approximation ratio γ is at most ∆+1· min{∆2 +2
1, |V |}. There are two possibilities to consider. In the ﬁrst case ∆ + 1 ≤ |V |. This
2
+1
implies ∆ = O( |V |) and γ ≤ ∆∆+1= O(∆) = O( |V |). In the second case|V |= O( |V |).
|V | < ∆2 + 1. This implies ∆ = Ω( |V |) and γ ≤ ∆+1
For many practical problems, such as meshes from ﬁnite element discretization ofPDEs, ∆2  |V |, making D2ColoringAlg an O(∆)-approximation algorithm.
The actual number of colors used by D2ColoringAlg depends on the order inwhich the vertices are visited. In our analysis, we assumed an arbitrary ordering.
A solution with fewer colors can be expected if a more elaborate ordering criterion
is used. For example, orderings such as largest ﬁrst, smallest last, incidence degree,
and saturation degree, appropriately adapted, are likely to be as fruitful for distance2 coloring as for distance-1 coloring. In the context of distance-1 coloring, theseorderings are analyzed in section 11.1. For a remark on the adaptation of theseorderings for distance-k coloring problems, see section 12.2.
3.5.2. A Partial Distance-2 Coloring Algorithm. D2ColoringAlg needs tobe modiﬁed only slightly to give an algorithm for the partial distance-2 coloring problem, our graph formulation of Problem 3.3. PartialD2ColoringAlg, given in
Algorithm 3.2, is a result of such a modiﬁcation.

652

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

Algorithm 3.2. A greedy partial distance-2 coloring algorithm.
procedure PartialD2ColoringAlg(Gb = (V1 , V2 , E))
Let v1 , v2 , . . . , v|V2 | be a given ordering of V2
Initialize forbiddenColors with some value a ∈ V2
for i ← 1 to |V2 | dofor each w ∈ N1 (vi ) dofor each colored vertex x ∈ N1 (w) do
forbiddenColors[color[x]] ← vi
end forend forcolor[vi ] ← min{c > 0 : forbiddenColors[c] = vi }end forend procedure

Let Gb = (V1 , V2 , E) be an input graph to this algorithm. For every vertexv ∈ V2 , the number of vertices at distance exactly two edges from v is at most∆(V2 ) · (∆(V1 ) − 1). Thus, vertex v can always be assigned a color from the set{1, 2, . . . , Cmax }, where Cmax = min{∆(V2 ) · (∆(V1 ) − 1) + 1, |V2 |}. Hence the arrayforbiddenColors in PartialD2ColoringAlg is of this size.
The time complexity result stated in Lemma 3.13 is evident.
Lemma 3.13. PartialD2ColoringAlg ﬁnds a partial distance-2 coloring of abipartite graph Gb = (V1 , V2 , E) on the vertex set V2 in time O(|V2 | · δ 1 (V2 ) · ∆(V1 )) =O(|E| · ∆(V1 )).
Let χ2 (Gb , V2 ) denote the minimum number of colors required for a partial distance2 coloring of Gb = (V1 , V2 , E) on V2 . One can then state the following analogous result
to Lemma 3.10.
Lemma 3.14. For every bipartite graph Gb = (V1 , V2 , E), χ2 (Gb , V2 ) ≥ ∆(V1 ).
In terms of Problem 3.3, Lemma 3.14 says that ∆(V1 ), which is equal to the
maximum number of nonzeros in a row of the Jacobian, is a lower bound on the least
number of groups required in a structurally orthogonal column partition. This lower
bound was observed by Coleman and Moré [30]. It should be noted that the size of alargest clique in Gb 2 [V2 ] is a tighter lower bound than ∆(V1 ).
3.5.3. A Distance-1 Coloring Algorithm. As mentioned earlier, Problem 3.3 canalso be formulated as a distance-1 coloring problem on the column intersection graph.
For comparison purposes, we consider D1ColoringAlg, the distance-1 analog of
D2ColoringAlg. The pseudocode for D1ColoringAlg (omitted here) diﬀers fromthat of D2ColoringAlg in only two ways: (1) the innermost for-loop is skipped;and (2) the array forbiddenColors is of size ∆ + 1, instead of min{∆2 + 1, |V |}.
The following result is straightforward.
Lemma 3.15. D1ColoringAlg ﬁnds a distance-1 coloring of a graph G =(V, E) in time O(|V |δ 1 ) = O(|E|).
From Lemmas 3.8, 3.13, and 3.15 it follows that the time required for constructing
and then distance-1 coloring the column intersection graph is asymptotically the sameas the time required for distance-2 coloring the column vertices of the bipartite graph.
This means that the two formulations are asymptotically comparable in terms ofoverall computation time. In practice, however, the runtimes while using the two
approaches may diﬀer considerably. Our experimental results reported in the next
subsection demonstrate this fact.

653

GRAPH COLORING FOR COMPUTING DERIVATIVES
Table 3.1 Matrix statistics.

Matrixlp cre a
lp dﬂ001
lp ken 11
lp stocfor3
lp ken 13
lp pds 10
lp maros r7
lhr10
lp pds 20
lp cre d
lp cre b
e30r2000
lhr14
lp ken 18
af23560
e40r0100
cage11lhr34
lhr71ccage12
lp osa 07
lp ﬁt2d

m
3,516
6,071
14,694
16,67528,632
16,5583,136
10,67233,874
8,926
9,648
9,661
14,270105,127
23,56017,281
39,08235,15270,354130,228
1,118
25

n
7,248
12,23021,34923,54142,65949,932
9,408
10,672
108,175
73,948
77,137
9,661
14,270154,699
23,56017,281
39,08235,15270,304130,228
25,067
10,524nnz
18,16835,632
49,05876,473
97,246
107,605
144,848
232,633
232,647
246,614260,785
306,356
307,858
358,171484,256
553,956
559,722764,014
1,528,092
2,032,536
144,812129,042

κmax
14
14
3
18
3
3
46
36
3
13
14
62
36
3
21
62
31
36
36
33
6
17

κmin
1
1
1
1
1
1
1
1
1
1
1
8
1
1
10
8
3
1
1
5
1
1

κavg
2.51
2.91
2.30
3.25
2.28
2.16
15.4
21.8
2.15
3.33
3.38
31.7
21.6
2.31
20.6
32.1
14.3
21.7
21.7
15.6
5.78
12.3

ρmax
360
228
122
15
170
96
48
63
96
808
844
62
63
325
21
62
31
63
63
33
17,613
10,500

ρmin
0
2
1
1
1
1
5
1
0
0
0
8
1
1
11
8
3
1
1
5
18
1,427

ρavg
5.17
5.87
3.34
4.59
3.40
6.50
46.2
21.8
6.87
27.6
27.0
31.7
21.6
3.40
20.6
32.1
14.3
21.7
21.7
15.6
130
5,162

3.6. Experimental Results. Here we report experimental results that compare
the bipartite graph–based distance-2 coloring formulation of Problem 3.3 with thecolumn intersection graph–based distance-1 coloring formulation. The comparison
focuses on the number of colors, overall execution time, and storage space requiredby the respective algorithms.
The algorithms in our test (including those reported in section 4.4) were implemented in the programming language C, and the experiments were conducted on anIntel Pentium 4, 2.53 GHz machine with 1 GB memory and 512 KB cache, runningLinux 2.4.20/RedHat 8.0.
The vertices are visited in their natural order when greedy coloring algorithms
are used in our experiments.
3.6.1. Test Matrices. Our testbed consists of matrices obtained from the University of Florida Sparse Matrix Collection [36]. Table 3.1 lists the test matrices alongwith some relevant structural statistics. Matrices with the preﬁx lp in their names arelinear programming problems translated from Netlib. (For brevity, we drop the preﬁx
lp when we refer to these matrices in later tables and ﬁgures.) The lhr -matrices arise
in steady-state chemical process simulation problems. The cage-matrices are models used in DNA electrophoresis. Matrices e30r2000 and e40r0100 come from ﬂuid
dynamics problems, and matrix af23560 is an airfoil problem.
In Table 3.1, the numbers of rows, columns, and nonzeros in each matrix are
listed under m, n, and nnz, respectively. The matrices are sorted in increasing orderof number of nonzeros. The maximum, minimum, and average number of nonzeros per column of a matrix are listed under κ with an appropriate subscript. The
corresponding ﬁgures per row are given under the various ρ’s.
3.6.2. Results and Discussion. Table 3.2 lists results obtained by running PartialD2ColoringAlg on the bipartite graph representation of the test matrices.

654

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN
Table 3.2 Results of partial d2 coloring the bipartite graph Gb (A).

Matrixcre adﬂ001
ken 11
stocfor3
ken 13
pds 10maros r7lhr10
pds 20cre d
cre b
e30r2000
lhr14
ken 18af23560
e40r0100
cage11lhr34
lhr71ccage12
Total1
osa 07ﬁt2d
Total2

|V |

|E|

∆

δ

δ

10,76418,301
36,04340,216
71,291
66,49012,544
21,344142,049
82,87486,78519,32228,540259,826
47,12034,562
78,16470,304
140,608
260,456
1,527,603
26,185
10,54936,734

18,16835,632
49,05876,473
97,246
107,605
144,848
232,633
232,647
246,614260,785
306,356
307,858
358,171484,256
553,956
559,722764,014
1,528,092
2,032,536
8,396,670
144,812129,042
273,854

360
228
122
18
170
96
48
63
96
808
844
62
63
325
21
62
31
63
63
33

0
1
1
1
1
1
1
1
0
0
0
8
1
1
10
8
3
1
1
5

3.38
3.89
2.72
3.80
2.73
3.23
23.1
21.8
3.28
5.95
6.01
31.7
21.6
2.76
20.6
32.1
14.3
21.7
21.7
15.6

17,613
10,500

1
1

11.1
24.5

K (ρmax )360 (360)
228 (228)130 (122)
16 (15)176 (170)
96 (96)74 (48)65 (63)96 (96)813 (808)
845 (844)
65 (62)65 (63)330 (325)
32 (21)66 (62)81 (31)65 (63)65 (63)96 (33)3,764 (3,573)
17,613 (17,613)
10,501 (10,500)
28,114 (28,113)

TGb

Tcol

Ttot

0
0
1
1
1
1
1
1
2
2
2
2
2
7
3
4
4
6
12
15
67
0
0
0

0.01
0.01
0.00
0.03
0.02
0.01
0.07
0.10
0.03
0.34
0.34
0.07
0.11
0.23
0.09
0.13
0.12
0.25
0.45
0.37
2.78
5.07
4.63
9.70

0.01
0.01
1.00
1.03
1.02
1.01
1.07
1.10
2.03
2.34
2.34
2.07
2.11
7.23
3.09
4.13
4.12
6.25
12.4
15.4
69.8
5.07
4.63
9.70

The left half of the table lists information on the underlying graph. The maximum,
minimum, and average degrees are given under ∆, δ, and δ, respectively. Note therelationship between the graph and matrix structural statistics given in Tables 3.1
and 3.2: |V | = m + n, |E| = nnz, ∆ = max{κmax , ρmax }, and δ = min{κmin , ρmin }.
The right half of Table 3.2 lists coloring and timing information. Each test matrix
is initially stored as an unordered list of nonzeros (edges). The time for reading
this data and for allocating memory for the various graph structures is not includedin the reported times. Column K in Table 3.2 lists the number of colors used byPartialD2ColoringAlg, and the number in parenthesis in the same column shows
the lower bound ρmax on the optimal number of colors. The time used (in seconds)
for constructing the bipartite graph (i.e., building the graph data structure from its
equivalent matrix data structure) is given in column TGb , and the time spent oncoloring is listed under Tcol . The last column gives the sum of the previous two. In
Table 3.2, as well as other later tables in this paper, if the time elapsed on somecomputation is less than 5 ms, we write the value 0.
Table 3.3 lists results obtained when D1ColoringAlg is run on the column
intersection graph of each test matrix. The column intersection graph is obtainedby ﬁrst constructing the bipartite graph and then computing the subgraph, inducedby the set of column vertices, of the square of the bipartite graph. See Lemma 3.7
for this relationship between the bipartite and column intersection graphs. In our
experiments, we used the routine ConstructG c, outlined in Algorithm 3.3, tocompute an intersection graph from a bipartite graph. In ConstructG c, the arraymarked, which is indexed with vertices in V2 , is used to avoid multiple edges between
two vertices in the intersection graph that are connected by several paths of length 2
in the bipartite graph.

655

GRAPH COLORING FOR COMPUTING DERIVATIVES
Table 3.3 Results of d1 coloring the column intersection graph Gc (A).

Matrixcre adﬂ001
ken 11
stocfor3
ken 13
pds 10maros r7lhr10
pds 20cre d
cre b
e30r2000
lhr14
ken 18af23560
e40r0100
cage11lhr34
lhr71ccage12
Total1
osa 07
ﬁt2d

|V |

|E|

∆

δ

δ

TGb

TGc

Tcol

Ttot

7,248
12,23021,34923,54142,65949,932
9,408
10,672
108,175
73,948
77,137
9,661
14,270154,699
23,56017,281
39,08235,15270,304130,228
930,536>>

253,411
250,976
459,921125,969
1,158,664
594,681
610,760
431,411
1,325,891
21,347,885
20,852,569
688,848
572,463
8,412,174
1,210,004
1,254,328
1,887,384
1,417,888
2,835,968
7,550,823
73,242,018
100
100

454
423
138
39
186
106
314
101
115
1,124
1,168
209
101
340
107
209
340
101
101
400

1
2
1
1
2
1
4
1
1
2
2
42
1
1
41
42
7
1
1
12

69.9
41.0
43.1
10.7
54.3
23.8
129
80.8
24.5
577
540
143
80.2
109
103
145
96.6
80.7
80.7
116

0
0
0
0
1
1
1
1
2
2
2
2
1
8
2
3
4
5
8
15
58

0
4
2
0
9
8
8
5
10
76
268
12
10
43
12
14
16
22
39
57
615

0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.02
0.02
0.25
0.75
0.02
0.00
0.15
0.02
0.01
0.02
0.03
0.04
0.11
1.51

0.01
4.01
2.01
0.01
10.0
9.01
9.01
6.02
12.0
78.2
271
14.0
11.0
51.1
14.0
17.0
20.0
27.0
47.0
72.1
675

millionmillion

edges
edges

Algorithm 3.3. An algorithm for computing the column intersection graph from abipartite graph.
procedure ConstructG c(Gb = (V1 , V2 , E))E ← ∅
for each v ∈ V2 dofor each w ∈ N1 (v) dofor each x ∈ N1 (w) do
if x = v thenif marked[x] = v thenE  ← E  ∪ (v, x)
marked[x] ← v
end ifend ifend forend forend forreturn Gc = (V2 , E  )end procedure

In Table 3.3 the time spent on the respective graph constructions is given undercolumns TGb and TGc . The other columns are deﬁned in a similar manner to Table 3.2.
Discussion. A comparison between the results in Tables 3.2 and 3.3 reveals,
among others, the following points.

656

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

Quality. The two coloring formulations are equivalent in terms of quality. In particular, as long as the vertices are visited in the same order, a greedy distance-1 coloringon a column intersection graph is exactly the same as a partial distance-2 coloringon a bipartite graph. In Table 3.3, the column listing the number of colors used by
the algorithm is omitted since the corresponding numbers are identical to those listed
under column K in Table 3.2. Another observation regarding the quality of coloringis that, as can be seen from Table 3.2, for the matrices in our experiment, the numberof colors used by the greedy algorithm is often close to, and in some cases even equal
to, the lower bound ρmax on the optimal value, hence implying optimality.
Graph Size. In general, for the matrices used in our experiments, a column intersection graph is larger in size than a bipartite graph and hence requires more storage
space. The diﬀerence in graph size becomes particularly large for matrices with relatively “dense” rows (“high” ρavg ). Using the quantity |V | + |E| as a measure of graphsize, for instance, the column intersection graph of matrix cre b is sixty times as large
as the corresponding bipartite graph. In the case of matrices osa 07 and ﬁt2d, thecolumn intersection graph was too big to ﬁt in the available memory space (the graphconstruction process in Algorithm 3.3 was stopped when the number of edges |E  |exceeded 100 million). Considering the sum total of the graph sizes of all other matrices from our testbed, the column intersection graph representation would require
more than seven times as much storage space as the bipartite graph representation.
Figure 3.2 shows a comparative plot of graph sizes using the two approaches for asubset of the computational results listed in Tables 3.2 and 3.3.
Time. The average time spent in constructing the bipartite graphs in our experiment is about 8 µs per edge, while the average time spent in constructing the column
intersection graphs is about 8.4 µs per edge. The partial distance-2 coloring algorithm
traverses on average |E|ρavg edges, while the distance-1 coloring algorithm traverses|E| edges, where |E| is the number of edges in the appropriate graph in each case.
The average time spent on coloring per edge is again roughly the same for the two algorithms with 17 ns for the distance-2 coloring algorithm and 20 ns for the distance-1
coloring algorithm.
In terms of overall computation time for each test case, our experimental results
show that a distance-2 coloring approach is signiﬁcantly faster than a method basedon distance-1 coloring (a subset of the results listed in Tables 3.2 and 3.3 is illustrated
in Figure 3.3). Again, the diﬀerence in overall execution time is large for matrices with
dense rows. For example, the overall time used for intersection graph construction
and distance-1 coloring for matrix cre b is 115 times that used for bipartite graphconstruction followed by partial distance-2 coloring.
For the ﬁrst twenty matrices, on which both approaches could be tested successfully, the total time spent on constructing and distance-1 coloring the column
intersection graph is about ten times the overall time spent on constructing anddistance-2 coloring the bipartite graph. This rather big diﬀerence is most likely due
to the fact that a column intersection graph–based approach requires larger memory
size, which in turn implies larger memory access time due to the multilevel memory
hierarchies in modern microprocessors [39]. A similar explanation can be oﬀered forthe observation that building a bipartite graph takes more time than distance-2 coloring it afterwards (see columns TGb and Tcol of Table 3.2). It should be pointed outthat in our experiments both the bipartite graph– and intersection graph–based approaches rely on “straightforward” implementations; the code was not been optimized
by performance-tuning.

657

GRAPH COLORING FOR COMPUTING DERIVATIVES

220|V | + |E|

214.22209.30

200
180
160
140
120
100
85.67
76.81

80
60
40
14.34
20
3.75
3.29
0 pds 20
cre d

29.1
22.9
19.26
15.3
14.5
12.34 12.71
8.3
6.38
5.31
5.89

3.47

6.18

cre b

ken 18 af23560 e40r0100 cage11

Bipartite graph

lhr34

lhr71ccage12MatrixColumn intersection graphFig. 3.2 Graph size comparison. The vertical axis shows |V | + |E| scaled by 100,000.

Time in sec.
315
270.75
300
285
270
255
240
225
210
195
180
165
150
135
120
105
78.25
72.11
90
75
51.15
47.04
60
45
20.02 27.03
17.01
30
12.45 15.37
7.23
6.25
4.13
4.12
2.34
15 2.340
cre d
cre b ken 18 e40r0100 cage11 lhr34 lhr71c cage12
Matrixd2-coloring Gb (A)d1-coloring Gc (A)Fig. 3.3 Overall execution time (graph construction plus coloring) comparison.

658

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

In practice, the time spent on distance-2 coloring a Jacobian is likely to be only asmall percentage of the time spent on computing the numerical values of the entriesof the Jacobian. For instance, a ﬁnite diﬀerence evaluation of a Jacobian matrix witha sparsity structure and size similar to that of matrix cage 12 takes about 4 seconds,
whereas coloring it takes less than 0.4 seconds [76].
4. Direct Computation of the Hessian. Given a twice continuously diﬀerentiable function f : Rn → R, the Hessian of f at the point x is the n × n symmetric∂ 2 f (x)matrix whose (i, j) entry ∇2 f (x)ij = ∂x
. When ∇f is available, ∇2 f can bei ∂xj
approximated by applying the ﬁnite diﬀerence formula (3.1) to the function F = ∇f .
We assume that evaluating the gradient ∇f (x) as a single entity is more desirable
than evaluating the components ∂1 f (x), . . . , ∂n f (x) separately. In this section, we
consider the eﬃcient computation of a Hessian matrix using a direct method.
4.1. The Matrix Partitioning Problem. Let di be a linear combination of asubset of the n coordinate vectors ej ∈ Rn . Analogous to the Jacobian estimationcase of the previous section, the problem of interest in the eﬃcient estimation of the
Hessian via a direct method can be stated in its general form as follows.
Problem 4.1. Given the sparsity structure of a symmetric matrix A, ﬁnd thefewest binary vectors d1 , d2 , . . . , dp such that the products Ad1 , Ad2 , . . . , Adp enable a
direct determination of A.
Powell and Toint [107] showed that in the case of Hessian estimation using ﬁnite
diﬀerences, in addition to exploiting sparsity using structurally orthogonal partitions,the number of function evaluations can be reduced further by exploiting symmetry.
Coleman and Moré [31] later deﬁned the following more general version of the partition
underlying the approach of Powell and Toint.
Definition 4.2. A partition of the columns of a symmetric matrix A is symmetrically orthogonal if for every nonzero element aij , either (1) the group containing thecolumn aj has no other column with a nonzero in row ri , or (2) the group containing
the column ai has no other column with a nonzero in row rj .
A more precise name for the partition referred to in Deﬁnition 4.2 would besymmetrically structurally orthogonal; we omit the word structurally for the sake ofbrevity.
The left part of Figure 4.1 illustrates a symmetrically orthogonal partition of thecolumns of a symmetric matrix. The ﬁgure shows a partition of the six columns of
the matrix into the three groups {a1 , a3 , a5 , a6 }, {a2 }, and {a4 }.
Let {C1 , C2 , . . . , Cp } be a symmetrically orthogonal partition of the columns ofa symmetric matrix A. With each group Ck , associate a binary vector dk having
components δj = 1 if aj belongs to Ck , and δj = 0 otherwise. Then
Adk =

n

j=1

δj aj =

aj .

aj ∈Ck

If aij = 0 and column aj is the only column in group Ck with a nonzero in row ri ,then aij = (Adk )i ; alternatively, if ai is the only column in a group Ck with a nonzero
in row rj , then aji = (Adk )j . This way, all the diagonal entries and at least one ofeach pair of symmetric entries of the matrix A can be determined directly with the pevaluations of Ad1 , . . . , Adp .
While using a ﬁnite diﬀerence technique to estimate a symmetric matrix A, if
a structurally orthogonal partition rather than a symmetrically orthogonal one isGRAPH COLORING FOR COMPUTING DERIVATIVES

1 2 3 4 5 6
X X
X X X

a1

659

a5
a2

X X

X X X
X

X X
X
X

Fig. 4.1

X
X

X

a6

a3
a4

A symmetrically orthogonal column partition and its representation as a star coloring of
the adjacency graph.

used, the estimate for aij may diﬀer from that of aji due to truncation error. Thus,
using a symmetrically orthogonal partition to estimate half of the oﬀ-diagonal nonzero
elements of a matrix and determining the other half by symmetry is preferable interms of both reducing computational work and ensuring that the computed matrix
is indeed symmetric.
Using Deﬁnition 4.2, Problem 4.1 can be cast as a partitioning problem as follows.
Problem 4.3. Given the sparsity structure of a symmetric matrix, ﬁnd a symmetrically orthogonal partition of its columns that has the fewest groups.
4.2. A Graph Coloring Formulation. The adjacency graph of a symmetric matrix exploits symmetry and hence is an appropriate graph representation for the sparsity structure of a Hessian. A Hessian matrix is often positive deﬁnite; hence theassumption in the deﬁnition of the associated adjacency graph that all diagonal entries of the matrix are nonzero is reasonable.
We proceed by presenting an adjacency graph–based characterization of structural
orthogonality in a symmetric matrix, a result due to McCormick [97].
Lemma 4.4. Let A 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 at a distance greater than two
from each other.
Proof. Assume that vertices ai and aj are at a distance greater than two fromeach other in G(A). This implies that (1) (ai , aj ) ∈ E, and (2) there exists no pathai , ak , aj in G(A) for any vertex ak ∈ V , 1 ≤ k ≤ n, k = i, j. In terms of matrix
A, this means that (1) aij = aji = 0, and (2) there is no index k ∈ [1, n], k = i, j,such that both aki and akj are nonzero. Recall that the diagonal entries aii and ajj
are nonzero. The ﬁrst statement above ensures that neither row index i nor j causes
columns ai and aj to be structurally nonorthogonal. The second statement says that
there is no other row index k that causes columns ai and aj to be nonorthogonal
either. Thus, columns ai and aj are structurally orthogonal.
Conversely, assume that columns ai and aj in A are structurally orthogonal.
Noting that aii = 0 and ajj = 0, this implies that (1) aij = aji = 0, and (2) there isno index k, 1 ≤ k ≤ n, k = i, j, such that aki = 0 and akj = 0. This means that inG(A) = (V, E), (1) (ai , aj ) ∈ E, and (2) there exists no path ai , ak , aj for any vertex
ak ∈ V , 1 ≤ k ≤ n, k = i, j. Hence, vertices ai and aj are at a distance d > 2 fromeach other. In particular, they are a distance d ≥ 3 apart from each other.
Notice that just as structural orthogonality of two columns in a nonsymmetricmatrix corresponds to the associated vertices in the bipartite graph being at a dis-

660

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

tance greater than two from each other, structural orthogonality of two columns in asymmetric matrix corresponds to the associated vertices in the adjacency graph being
at a distance greater than two from each other. There is, however, one diﬀerence
between the two cases: the distance between a pair of column vertices in a bipartite
graph is always even, while the distance between a pair of vertices in an adjacency
graph could be odd or even.
By Lemma 4.4, ﬁnding a structurally orthogonal partition of the columns of asymmetric matrix A is equivalent to ﬁnding a distance-2 coloring of the adjacency
graph G(A).
We note that Lemma 4.4 is also applicable to a Jacobian matrix that is symmetricin its zero-nonzero structure, but not in the numerical values of its entries. In particular, distance-2 coloring the adjacency graph of a structurally symmetric Jacobianmatrix with nonzero diagonal entries is equivalent to a structurally orthogonal column
partition.
In the case of Hessian estimation, the fact that entries aij and aji of matrix A
are equal in numerical value can be exploited to further reduce the number of groups(colors) required. We now consider the graph coloring formulation of Problem 4.3,
the partitioning problem where symmetry is exploited.
Consider a symmetric matrix A with nonzero diagonal elements and let aij , i = j,be any nonzero element in A. Further, let aki , k = i, j, and ajl , l = i, j, k, be anyother two nonzero elements. These three nonzeros (and their symmetric counterparts)and the relevant diagonal entries of the matrix A have the following structure:


aii
 aji

 aki

aij
ajj

aik
akk

aj


aj 
.

a

By Deﬁnition 4.2, in a symmetrically orthogonal partition of the columns of A,
R1. each of the column pairs (ai , aj ), (ai , ak ), and (aj , al ) have their respective
elements in two diﬀerent groups; andR2. either columns aj and ak belong to two diﬀerent groups, or columns ai andal belong to two diﬀerent groups.
In fact, these two requirements are also suﬃcient for a symmetrically orthogonal
partition. Notice that a partition of the four columns ai , aj , ak , and al into the threegroups {ai }, {aj , ak }, and {al } satisﬁes R1 and R2.
Coleman and Moré [31] characterized the requirements R1 and R2 in terms of acoloring in the associated adjacency graph.
Definition 4.5. A mapping φ : V → {1, 2, . . . , p} is a star coloring of the graph
G = (V, E) if (1) φ is a distance-1 coloring of G, and (2) every path on four vertices
uses at least three colors.
Notice that R1 corresponds to the ﬁrst requirement of Deﬁnition 4.5, and R2
corresponds to the second.
The subﬁgure on the right in Figure 4.1 shows a star coloring representation of
the symmetrically orthogonal partition shown on the left. In Figure 4.1, vertices a1 ,
a3 , a5 , and a6 are assigned one color; vertex a2 is assigned a second color; and vertex
a4 is assigned a third color.
Coleman and Moré [31] refer to the coloring given in Deﬁnition 4.5 as path coloring. Unfortunately, the term path coloring is used in the optical networks literature torefer to an assignment of colors to paths in a graph serves as a theoretical model forGRAPH COLORING FOR COMPUTING DERIVATIVES

661

routing problems. Examples of works that discuss such applications include [1, 108],
and a recent survey is available in [15]. The term star coloring owes its name to thefact that in such a coloring, every subgraph induced by vertices assigned any two colors is a collection of stars. As an illustration, consider the graph shown in Figure 4.1
and the subgraph induced by vertices a1 , a2 , a3 , a5 , and a6 therein.
The following theorem formalizes the connection between symmetrically orthogonal partitioning and star coloring. The result follows from the discussion that led tothe deﬁnition of star coloring.
Theorem 4.6 (Coleman and Moré [31]). Let A be a symmetric matrix with
nonzero diagonal elements, and let G(A) = (V, E) be its adjacency graph representation. A mapping φ is a star coloring of G(A) if and only if φ induces a symmetricallyorthogonal partition of the columns of A.
By Theorem 4.6, the following problem is equivalent to Problem 4.3.
Problem 4.7. Given the adjacency graph G(A) = (V, E) representing the sparsity structure of a symmetric matrix A having nonzero diagonal elements, ﬁnd a starcoloring of G(A) that uses the fewest colors.
4.3. Algorithms. Coleman and Moré [31] showed that the problem of ﬁnding astar coloring with the fewest colors is NP-hard even if the graph is bipartite. Here we
discuss two heuristic algorithms for this problem. The ﬁrst algorithm is new, while thesecond was proposed earlier by Powell and Toint [107] in terms of matrices. In ﬁnding
a valid color to assign a vertex, the ﬁrst algorithm visits the distance-3 neighbors ofthe vertex while the second algorithm visits only the distance-2 neighbors. The latteressentially solves a more restricted coloring problem. The two algorithms represent
a trade-oﬀ between the number of colors used and runtime. In both algorithms anarray forbiddenColors of size Cmax = min{∆2 + 1, |V |} is used.
While this paper was in print, we developed a new star coloring algorithm whosetime complexity is the same as the second algorithm discussed here and whose performance is similar to the ﬁrst. See section 12.2 for a remark.
4.3.1. The First Star Coloring Algorithm. StarColoringAlg1 (Algorithm
4.1) outlines the ﬁrst algorithm. Figure 4.2 graphically shows the decision madeduring one of the |V | steps of the for-loop in lines 4–21 of StarColoringAlg1. In
Figure 4.2, 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 tree nodes at levels w, x, and y, respectively. A green-painted node signiﬁesthat 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
y. The correspondence between Figure 4.2 and lines 6, 10, and 14 of Algorithm 4.1 isobvious. Note that each tree node corresponds to many vertices of the input graph.
Line 6 in Algorithm 4.1 guarantees that the resulting coloring is consistent with adistance-1 coloring. Notice that in line 10 of the algorithm, the color of vertex x in apath v, w, x where w is not yet colored is forbidden for vertex v. Later on, when vertex
w is colored, the test in line 6 ensures that w will not have two neighbors of the samecolor and thus will not be the center of a 2-colored path of length two (which could
stretch into a 2-colored path of length three). Thus it remains to ensure that v does not
become an endpoint of a 2-colored path of length three, and this is done in line 14 of thealgorithm. The break statement in Algorithm 4.1 is used since the discovery of the ﬁrst
(out of possibly several) 2-colored paths of length two connected to v suﬃces to force
v to have a third color. (A break statement forces control to jump out of the nearest
loop that contains the statement, in this case, out of the for-loop in lines 13–18.)

662

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

Algorithm 4.1. A greedy star coloring algorithm.
1: procedure StarColoringAlg1(G = (V, E))
Let v1 , v2 , . . . , v|V | be a given ordering of V
2:
3:
Initialize forbiddenColors with some value a ∈ V
4:for i ← 1 to |V | do5:
for each w ∈ N1 (vi ) do
6:
if w is colored thenforbiddenColors[color[w]] ← vi
7:
8:
end if9:
for each colored vertex x ∈ N1 (w) do
10:
if w is not colored thenforbiddenColors[color[x]] ← vi
11:
12:
else
13:
for each colored vertex y ∈ N1 (x), y = w do
14:
if color[y] = color[w] thenforbiddenColors[color[x]] ← vi
15:
16:
break
end if
17:
18:
end forend if
19:
20:
end forend for21:
22:
color[vi ] ← min{c > 0 : forbiddenColors[c] = vi }end for23:
24: end procedure

v

f
w
f

?
x

y

Fig. 4.2 Visualizing a step in StarColoringAlg1.

The work done to assign a color to vertex vi in StarColoringAlg1 is proportional to d3 (vi ). Thus we get the following result.
Lemma 4.8. StarColoringAlg1 ﬁnds a star coloring of a graph G = (V, E)
in time O(|V |δ 3 ).
4.3.2. The Second Star Coloring Algorithm. Recall that a coloring of a graphis an assignment of positive integers to its vertices. The second star coloring algorithm
is based on the observation that a star coloring is a relaxed distance-2 coloring, andmakes use of the fact that colors are positive integers.
One way of relaxing the requirement for a distance-2 coloring so as to obtain astar coloring is to let two vertices at distance of exactly two edges from each otherGRAPH COLORING FOR COMPUTING DERIVATIVES

663

Algorithm 4.2. A second greedy star coloring algorithm.
1: procedure StarColoringAlg2(G = (V, E))
Let v1 , v2 , . . . , v|V | be an ordering of V
2:
3:
Initialize forbiddenColors with some value a ∈ V
4:for i ← 1 to |V | do5:
for each w ∈ N1 (vi ) do
6:
if w is colored thenforbiddenColors[color[w]] ← vi
7:
8:
end if9:
for each colored vertex x ∈ N1 (w) do
10:
if w is not colored thenforbiddenColors[color[x]] ← vi
11:
12:
else
13:
if color[x] < color[w] then14:
forbiddenColors[color[x]] ← vi
end if15:
16:
end ifend for17:
18:
end forcolor[vi ] ← min{c > 0 : forbiddenColors[c] = vi }19:
20:
end for
21: end procedure

share a color as long as the vertex in between them has a color of lower value. Moreprecisely, let v, w, x be a path in G and suppose v and w are colored and we want
to determine the color of x. Clearly, we need to make sure that φ(x) is distinct
from φ(w). Further, we allow φ(x) to be equal to φ(v) as long as φ(w) < φ(v). To
see that this coloring can always be extended to yield a valid star coloring, consider
the path v, w, x, y, an extension of the path v, w, x in one direction. Now, sinceφ(x) = φ(v) > φ(w), we cannot let φ(y) be equal to φ(w). Obviously, φ(y) should
be diﬀerent from φ(x), otherwise it will not be a valid distance-1 coloring. Thus thepath v, w, x, y uses three colors, φ is a distance-1 coloring, and therefore it is a valid
star coloring. StarColoringAlg2, the algorithm that uses this idea, is outlined in
Algorithm 4.2. The work done to assign a color to vertex vi in StarColoringAlg2
is proportional to d2 (vi ), hence the following result.
Lemma 4.9. StarColoringAlg2 ﬁnds a star coloring of a graph G = (V, E)
in time O(|V |δ 2 ).
Notice that the coloring produced by StarColoringAlg2 is a more restrictedvariant of star coloring, and therefore StarColoringAlg1 is likely to use fewer colors than StarColoringAlg2 for a given input graph. Figure 4.3 shows an example
where the ﬁrst algorithm uses three colors in coloring the vertices in their alphabeticalorder while the second uses four in doing the same. The two star coloring algorithms
represent a trade-oﬀ between the number of colors used and the computational timeneeded.
4.4. Experimental Results. Our experimental work in the Hessian estimationcase has two objectives: (i) to experimentally demonstrate the advantage of exploiting
symmetry, and (ii) to compare and contrast the performance of the two star coloringalgorithms.

664

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

b 2

Fig. 4.3

b 2

1

3

1

1

3

4

a

c

d

a

c

d

A simple example where StarColoringAlg1 uses fewer colors than StarColoringAlg2.
Table 4.1 Graph statistics.

Graph
mrng1
mrng2
598a
144
m14b
auto

Table 4.2

|V |257,000
1,017,253
110,971
144,649
214,765448,695|E|505,048
2,015,714
741,934
1,074,393
1,679,018
3,314,611

∆
4
4
26
26
40
37

δ
2
2
5
4
4
4

δ
3
3
13
14
15
14

The number of colors used by and the runtimes of D2ColoringAlg, StarColoring
Alg1, and StarColoringAlg2.

Graph
mrng1
mrng2
598a
144
m14b
auto
Total

K(d2)
12
12
38
41
42
42
187

K(star1)
8
9
27
28
29
29
130

K(star2)
10
10
32
35
34
36
157

T (d2)
0.37
1.74
0.74
0.90
1.19
3.97
8.91

T (star1)
0.61
2.90
3.00
4.35
5.65
17.4
33.9

T (star2)
0.35
1.75
0.75
0.93
1.23
4.07
9.08

4.4.1. Test Graphs. The matrices behind the test graphs used in our experiments
arise from ﬁnite element methods [46]. Table 4.1 gives some structural information
about the adjacency graphs of these matrices.
4.4.2. Results and Discussion. Table 4.2 shows the performance of algorithms
D2ColoringAlg, StarColoringAlg1, and StarColoringAlg2 when applied
to our test graphs. In each algorithm, the vertices are visited in their natural order.
The left half of the table shows the number of colors used by the diﬀerent algorithms
and the right half shows the corresponding time (in seconds) spent on coloring.
The results clearly demonstrate the advantage of exploiting symmetry: star coloring requires signiﬁcantly fewer colors than distance-2 coloring. The total numberof colors required by StarColoringAlg1 over all involved test graphs is about
30% less than the number required by D2ColoringAlg. The table also shows
nicely the time/quality trade-oﬀ between the two star coloring algorithms. StarColoringAlg2 uses nearly the same time as D2ColoringAlg, while the numberof colors required is 16% less.
5. Bidirectional, Direct Computation of the Jacobian. Recall that a bidirectional partition, as opposed to a unidirectional partition, involves both the rows andcolumns of a matrix. We begin this section by motivating the need for a bidirectional

GRAPH COLORING FOR COMPUTING DERIVATIVES

665

partition within the context of AD. We then introduce the partitioning problem, discuss its graph coloring formulation, and suggest algorithms for solving the problem.
5.1. The Need for a Bidirectional Partition. Through its forward and reverse
mode AD allows a columnwise and a rowwise computation of a Jacobian, respectively.
One way of exploiting sparsity in a matrix in this context is to separately partition thecolumns and rows of the matrix, choose the partition which gives the minimum numberof groups, and apply the appropriate mode of AD to compute the matrix entries. Fora symmetric matrix, a row partition is equivalent to a column partition, but fora nonsymmetric matrix, the two partitions may diﬀer considerably. For example,
consider an n × n matrix where all the entries on the diagonal and the ﬁrst row arenonzero, and the rest of the matrix entries are all zero. In such a case, a structurally
orthogonal column partition requires n groups, whereas a structurally orthogonal 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 arrowhead matrix where all of the elementsin the ﬁrst row, ﬁrst column, and the diagonal are nonzero and the rest of the entries
are all zero. For such a pattern, a structurally orthogonal row partition requires ngroups and so does a column partition. However, using a combined row and column
partition, three groups are enough to determine all the nonzero entries of the matrix.
First, separately evaluate the entries in the ﬁrst column and the ﬁrst row (two groups).
Then, since the remaining (n − 1) × (n − 1) matrix is diagonal, group the columns
together and determine all entries by one forward AD pass. Thus, three groups (twocolumn and one row) suﬃce to determine all the nonzero entries.
In the current section we consider such a computation of a nonsymmetric matrix
using the combined modes of AD via a direct method. The corresponding problem
using a substitution method will be discussed in section 6.2.
Note that a bidirectional partition does not make sense for computing a numerically symmetric matrix. In particular, a symmetry-exploiting unidirectional partition,
as discussed in section 4, is suﬃcient.
5.2. The Matrix Partitioning Problem. Let A be an m × n matrix with a known
sparsity structure. Partitioning the columns of A into p groups can be seen as thetask of seeking an n × p binary matrix D whose (j, k) entry is deﬁned as follows:
(5.1)

djk =

1 if column aj belongs to group k,
0 otherwise.

Recall that D has been referred to as a seed matrix in earlier sections. Using this,an alternative way of posing Problem 3.1 would then be the following: Given the
structure of an m × n matrix A, ﬁnd an n × p seed matrix D with the least value ofp such that the product AD enables a direct determination of A. By the same token,
the problem that arises in the bidirectional computation of a Jacobian via a direct
method can be posed as follows.
Problem 5.1. Given the sparsity structure of an m × n matrix A, ﬁnd an
n × p1 binary matrix D1 and an m × p2 binary matrix D2 such that the products AD1
and D2T A together enable a direct determination of A and the value p = p1 + p2 is
minimized.
In this problem formulation, the product AD1 corresponds to the columns of A
determined using the forward mode of AD, and the product D2T A corresponds to therows of A determined using the reverse mode of AD.

666

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

1
2
3
4
5

X

X

X X
X

X
X

X X X X
X

X

1 2 3 4 5
Fig. 5.1

A structurally orthogonal bidirectional partition of a matrix. The color(s) used on eachnonzero entry shows the column and/or row group from which it can be computed.

Hossain and Steihaug [62] studied Problem 5.1 and identiﬁed an associated partitioning problem. Their formulation relies on the notion of a consistent row-column
partition in which the entire set of rows and columns is partitioned into two respective
sets of groups. Coleman and Verma [32] also studied the same problem and identiﬁed
a similar bidirectional partitioning problem. Their notion of partition diﬀers fromthat of Hossain and Steihaug in that it partitions only a subset of the columns andthe rows of the matrix that suﬃces for the direct determination of the entries. We
consider such a partition and introduce its formal deﬁnition below.
Definition 5.2. A bidirectional partition of a matrix A is a pair (ΠC , ΠR ) whereΠC is a partition of a subset of the columns of A and ΠR is a partition of a subset ofthe rows of A.
Definition 5.3. A bidirectional partition (ΠC , ΠR ) of a matrix A is structurally
orthogonal if for every nonzero element aij , either (1) column aj is in a group of ΠC
which has no other column having a nonzero in row ri , or (2) row ri is in a group of
ΠR which has no other row having a nonzero in column aj .
The number of column and row groups in a structurally orthogonal bidirectional
partition corresponds to the number of forward and reverse AD passes, respectively,
required to compute the nonzero entries directly. To see this, observe that a nonzero
element aij can be determined either from a column group where column aj is the
only column with a nonzero in row ri , or from a row group where row ri is the only
row with a nonzero in column aj .
Figure 5.1 shows an example of a structurally orthogonal bidirectional partitionof a matrix. In the example, ΠC includes columns a1 and a3 , whereas ΠR includes
rows r2 , r4 , and r5 . As can be seen from the color used at the left and bottom
edges of the ﬁgure, column a1 forms one column group and column a3 forms another;similarly, rows r2 and r5 form one group and row r4 forms a second row group. Thus,
the bidirectional partition uses a total of four groups. Notice that some entries of thematrix (painted with one color) can be computed only from one group while others(painted with two colors) can be computed from either of the two groups. As thereader can easily verify, for this example, a row-only or a column-only structurally
orthogonal partition would have required ﬁve groups.
Assuming that the computational costs involved in the forward and reverse modesof AD are of the same order, in an eﬃcient method that uses a bidirectional partition (ΠC , ΠR ), the value |ΠC | + |ΠR | is required to be as small as possible. Thus
Problem 5.1 can be cast as a partitioning problem in the following way.

667

GRAPH COLORING FOR COMPUTING DERIVATIVES

1
2
3
4
5

X

X

X X
X

X
X

X X X X
X

X

r1
r2
r3
r4
r5

a1
a2
a3
a4
a5

1 2 3 4 5

Fig. 5.2

A structurally orthogonal bidirectional partition and its representation as a star bicoloring.
The illustration shows a bidirectional partition (ΠR , ΠC ): ΠR = {r2 , r5 }; {r4 } and ΠC ={a1 }; {a3 }.

Problem 5.4. Given the sparsity structure of an m × n matrix A, ﬁnd a structurally orthogonal bidirectional partition (ΠC , ΠR ) such that |ΠC |+|ΠR | is minimized.
It should be pointed out that in formulating Problem 5.4 within the context ofAD, we are concerned only with computational cost. However, in general, the forward
mode of AD requires less memory space than the reverse mode, making the former
perhaps more desirable. Hence, a more accurate objective would be to minimize
w1 |ΠC | + w2 |ΠR | for some empirically determined weights w1 and w2 [66].
5.3. A Graph Coloring Formulation. When a nonsymmetric matrix A is represented by its bipartite graph Gb (A) = (V1 , V2 , E), we have seen that a structurally
orthogonal unidirectional partition can be obtained by ﬁnding a partial distance-2
coloring of Gb on the set V2 of column vertices. We now consider how this coloringneeds to be modiﬁed to ﬁnd a structurally orthogonal bidirectional partition. Noticethat the coloring we are looking for has to meet the following conditions.
• Some vertices may not be involved in the determination of any nonzero entryof the underlying matrix. Such vertices are assigned the “neutral” color zero.
We use positive integers to denote the other colors.
• The colors assigned to vertices in V1 should be disjoint from colors assigned
to vertices in V2 , except for the neutral 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 positive color.
• A nonzero matrix entry may be determined from either a positively coloredcolumn vertex or a positively colored row vertex. This suggests that thecoloring condition sought here is some relaxation of the distance-2 coloring
requirement imposed in the case of unidirectional partition.
The following deﬁnition, introduced by Coleman and Verma [32], makes the aforementioned conditions more precise. The subsequent theorem establishes the equivalence between the matrix and graph problems.
Definition 5.5. Let Gb = (V1 , V2 , E) be a bipartite graph. A mapping φ :[V1 , V2 ] → {0, 1, . . . , p} is a star bicoloring of Gb if the following conditions hold:
1. If u ∈ V1 and v ∈ V2 , then φ(u) = φ(v) or φ(u) = φ(v) = 0.
2. If (u, v) ∈ E, then φ(u) = 0 or φ(v) = 0.
3. If vertices u and v are adjacent to a vertex w with φ(w) = 0, then φ(u) = φ(v).
4. Every path on four vertices uses at least three colors.
Figure 5.2 shows a structurally orthogonal bidirectional partition of a matrix andits representation as a star bicoloring in the associated bipartite graph. The reader

668

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

Algorithm 5.1. A scheme for star bicoloring.
procedure StarBicoloringScheme(Gb = (V1 , V2 , E))
1. Find a suitable vertex cover C in Gb
2. Assign the vertices in the set I = (V1 ∪ V2 ) \ C the color 0
3. Color the vertices in C such that the result is a star bicoloring of Gb
end procedure

could verify that these four conditions listed above are satisﬁed by the coloring in thegraph.
Theorem 5.6 (Coleman and Verma [32]). 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 star
bicoloring of Gb (A) if and only if φ induces a structurally orthogonal bidirectional
partition (ΠC , ΠR ) of A.
Thus Problem 5.4 is equivalent to the following graph coloring problem.
Problem 5.7. Given the bipartite graph Gb (A) = (V1 , V2 , E) representing thesparsity structure of an m × n matrix A, ﬁnd a star bicoloring of Gb (A) that uses the
fewest colors.
5.4. Algorithms. In a star bicoloring of a bipartite graph, some vertices are
assigned the neutral color zero. (As an example, see the unpainted vertices in Figure 5.2.) We make the following observation which helps us identify a possible setof such vertices. The observation is a direct consequence of conditions 1 and 2 ofDeﬁnition 5.5.
Observation 5.8. Let Gb = (V1 , V2 , E) be a bipartite graph and φ : [V1 , V2 ] →{0, 1, . . . , p} be a star bicoloring of Gb . Then
• the set C = {v : φ(v) = 0} is a vertex cover in Gb , and• the set I = {v : φ(v) = 0} is an independent set in Gb .
One consequence of Observation 5.8 is that |I|+|C| = |V1 |+|V2 |. Thus, a decreasein the cardinality of the vertex cover C results in an increase in the cardinality of the
independent set I.
Observation 5.8 suggests a scheme for solving Problem 5.7. This scheme, called
StarBicoloringScheme, is outlined in Algorithm 5.1.
In step 1 of starBicoloringScheme, a vertex cover needs to be chosen carefully. Though a vertex cover of small size is desirable, minimizing its cardinality is
not a primary objective. Speciﬁcally, the chosen vertex cover should be such that it
results in fewer colors compared to both a partial distance-2 coloring on V1 (rowwise
unidirectional partition) and a partial distance-2 coloring on V2 (columnwise unidirectional partition). As the discussion in section 5.1 suggests, such a vertex coverneeds to include vertices from V1 and V2 having a relatively high number of distance1 neighbors, even when the implication is that the vertex cover is not minimal in size.
Coleman and Verma [32] have suggested a procedure (formulated in matrix terms)
for computing a vertex cover suitable for star bicoloring.
A high-level translation of the procedure of Coleman and Verma in graph terms
would be as follows. A vertex cover C = C1 ∪ C2 is obtained by ﬁrst ﬁnding an
independent set I = I1 ∪ I2 and then applying the diﬀerences C1 = V1 \ I1 andC2 = V2 \ I2 . The independent set is obtained using an algorithm that has the
ﬂavor of a minimum degree algorithm for ﬁnding a maximal independent set. An
additional ingredient of the algorithm is a classiﬁcation of the edge set E into four
categories. The ﬁrst category corresponds to edges not yet covered. An edge in the

669

GRAPH COLORING FOR COMPUTING DERIVATIVES

v

w
f

f

?
x

y

Fig. 5.3 Visualizing a step in StarBicoloringAlg.

second category is an edge covered solely by a vertex from V1 , an edge in the third
category is covered solely by a vertex from V2 , and an edge in the fourth category iscovered both by a vertex from V1 and by a vertex from V2 . This classiﬁcation is used
to deﬁne specialized vertex “degrees” (number of incident edges of each category).
Initially all edges are set to be in category 1, to reﬂect that they are not covered yet.
In each step of the greedy algorithm, a vertex v1 of minimum degree in V1 and avertex v2 of minimum degree in V2 are chosen; the minimum degree in each case isevaluated using the subgraph of Gb obtained by deleting category-2 and category-3
edges. Then, based on a test function that involves a comparison between current
maximum degrees (restricted to category-2 and category-3 edges), one of the vertices
v1 or v2 is chosen to be included in the independent set. Once the vertex to be includedin the independent set is determined, the edges incident on the chosen vertex and its
neighbors are placed into appropriate categories that reﬂect the choice. The algorithm
terminates when no category-1 edge is left (all edges are covered).
Once steps 1 and 2 are carried out, step 3 can be done by a suitable adaptation of
StarColoringAlg1. Let StarBicoloringAlg be such an adaptation. One of thediﬀerences between the coloring and bicoloring algorithms is that in the latter case,two disjoint sets of colors are used in coloring the vertices in V1 and V2 of the bipartite
graph Gb = (V1 , V2 , E). Another diﬀerence is that at the instance of the bicoloringalgorithm in which vertex v is colored, a vertex within the distance-3 neighborhoodof 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 vertex v thus needs to consider thesethree options.
Figure 5.3 shows a visual presentation of an iteration (in the loop over all vertices)of StarBicoloringAlg. Note the similarity with Figure 4.2. Since the colors for thevertices in V1 and V2 are chosen from two disjoint sets, in choosing a color for vertex v
in Gb = (V1 , V2 , E), we need only consider colors of vertices that are exactly two edges
away from v. In the ﬁgure, green-painted nodes correspond to vertices with positive
colors, shaded nodes show vertices with color zero, and unpainted nodes correspondto uncolored vertices. Observe that the node with color zero at level w has only twochildren; it cannot have a child with color zero, for otherwise there would exist an
uncovered edge. The colors of the vertices in the nodes marked by an f indicate forbidden colors; whether the color at the node marked by “?” is forbidden or not depends
on the color used at node y: if φ(w) = φ(y), φ(x) is forbidden; otherwise, it is not.
The time complexity of StarBicoloringAlg is O((|V1 | + |V2 |)δ 3 ), which is also
the overall time complexity of StarBicoloringScheme assuming that step 1 is doneusing a greedy algorithm that is linear in the number of edges.

670

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

Notice that a partial distance-2 coloring of Gb on V2 is just a special case of
StarBicoloringScheme. To see this, consider the trivial choice of vertex cover
C = V2 in step 1. This implies that, in step 2, the vertices in the set I = V1 will becolored with color zero. By condition 3 of Deﬁnition 5.5, vertices adjacent to a vertex
colored with color zero are required to be assigned diﬀerent colors. Thus, the result
is eﬀectively a partial distance-2 coloring of Gb on V2 .
Hossain and Steihaug [62] and Coleman and Verma [32] each proposed an algorithm for Problem 5.7. These can be interpreted in light of StarBicoloringScheme.
The algorithm of Hossain and Steihaug implicitly ﬁnds a vertex cover while the coloring of the graph proceeds. Using our terminology, the vertices that remain uncoloredat the end of the Hossain–Steihaug algorithm form an independent set in the graphand can thus be assigned the neutral color zero.
The algorithm of Coleman and Verma uses a preprocessing step to identify therows and columns of the underlying matrix that eventually need to be colored withpositive values. The preprocessing step uses a matrix-based procedure which eﬀectively produces a small-sized vertex cover. After the preprocessing step, certain “column and row intersection graphs,” adapted to the star bicoloring requirements, areconstructed to ﬁnally use known distance-1 coloring heuristics on the resulting graphs.
6. Substitution Methods. In a unidirectional computation of a matrix A via
a substitution method, the vectors d1 , d2 , . . . , dp are chosen such that the systemof equations deﬁned by the products Ad1 , Ad2 , . . . , Adp can be made triangular. In
this case, the partition deﬁned by the p vectors could fulﬁll a more relaxed set of
requirements compared to a structurally orthogonal partition used in a direct method,and hence results in a smaller p. This fact has been especially useful when estimatinga Hessian matrix since substitution can be eﬀectively combined with the exploitationof symmetry [27, 31]. In the Jacobian case, a substitution method is worth considering
when the computation is bidirectional [32].
In this section, within the context of substitution methods, we consider unidirectional partitions for Hessians (section 6.1) and bidirectional partitions for Jacobians(section 6.2).
An example of a variant of a substitution method for a unidirectional computationof a Jacobian has been considered by Hossain and Steihaug [63]. The method relies
on ﬁrst ﬁnding a structurally orthogonal partition of the columns and then merging apair of consecutive groups to get a column grouping that allows overlaps. In this way
a structurally orthogonal partition consisting of p groups is changed into a column
grouping having p − 1 groups. The former deﬁnes a diagonal system of equations,while the latter leads to a triangular system.
6.1. Computing the Hessian. To illustrate the fact that a partition used in a
substitution method requires fewer groups than one used in a direct method, consider
the 4 × 4 symmetric matrix A shown below:

a11
a21



a12
a22
a32


a23
a33
a43


.
a34 
a44

Any symmetrically orthogonal partition of the columns of this matrix—and hencea direct method—requires at least three groups. An example of a symmetricallyorthogonal partition is {a1 , a3 }, {a2 }, and {a4 }. However, if we do not insist on

GRAPH COLORING FOR COMPUTING DERIVATIVES

671

determining the elements directly, two groups would suﬃce. For example, consider
the partition {a1 , a3 } and {a2 , a4 }. The seed matrix corresponding to this partitionis
1 0 1 0.
DT =0 1 0 1The product AD from which the nonzeros of A can be recovered is


a11
a12
a12 + a23a22 .
AD = 
 a33
a23 + a34 a34
a44
Thus, nonzero entries a11 , a12 , a22 , a33 , a34 , and a44 can be obtained directly andelement a23 can be obtained via substitution either from the expression a12 + a23
or from a23 + a34 . Note that in the product AD, every nonzero entry aij has been
identiﬁed with its symmetric counterpart aji .
6.1.1. The Matrix Partitioning Problem. In general, a column partition of asymmetric matrix induces a substitution method if there is an ordering of the nonzero
entries (unknowns) such that all unknowns can be solved for, in that order, using
symmetry and previously solved elements.
We formally deﬁne such a partition and then state the corresponding partitioning
problem. In the following deﬁnition, we identify the two matrix elements aij andaji of the symmetric matrix A. The ordering of the elements in the deﬁnition is theordering in which the matrix elements are evaluated in a substitution method.
Definition 6.1. A partition of the columns of a symmetric matrix A is said tobe substitutable if there exists an ordering on the elements of A such that for every
nonzero aij , either (1) column aj is in a group where all the nonzeros in row ri , fromother columns in the same group, are ordered before aij , or (2) column ai is in a groupwhere all the nonzeros in row rj , from other columns in the same group, are ordered
before aij .
Problem 6.2. Given the sparsity structure of a symmetric matrix, ﬁnd a substitutable partition of its columns that has the fewest groups.
6.1.2. A Graph Coloring Formulation. To formulate Problem 6.2 as a graphproblem, we need to introduce the notion of acyclic coloring.
Definition 6.3. A mapping φ : V → {1, 2, . . . , p} is an acyclic coloring of agraph G = (V, E) if (1) φ is a distance-1 coloring, and (2) every cycle in G uses atleast three colors.
Note that since the vertices in every cycle are assigned at least three colors in an
acyclic coloring, the subgraph induced by the set of vertices assigned any two colors
does not contain a cycle, and hence it is a forest. The name acyclic coloring stems
from this fact.
Consider the 10 × 10 symmetric matrix A whose nonzero structure is depicted in
Figure 6.1. The ﬁgure shows a partition of the columns of A into the three groups{a1 , a3 , a5 }, {a2 , a4 , a7 , a9 }, and {a6 , a8 , a10 }. We will show that this partition issubstitutable.
An acyclic coloring of the adjacency graph G(A) is shown in Figure 6.1 (right);
the vertices {a1 , a3 , a5 } are colored red, the vertices {a2 , a4 , a7 , a9 } are colored blue,
and the vertices {a6 , a8 , a10 } are colored green. Note that the subgraph induced

672

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

a1

1 2 3 4 5 6 7 8 9 10
XX
X
XXX X
XXX X
XX
X
X
XX X
X XX
X
X
XX
X XXX
X XXX
X
XX

a3

a5

a7
Fig. 6.1

a2

a4

a6

a8

a9

a 10

A substitutable partition of the columns of a symmetric matrix and its representation as
an acyclic coloring of the associated adjacency graph.

by the vertices colored red or blue, the set {a1 , a2 , a3 , a4 , a5 , a7 , a9 }, is a forest; thevertex a9 is an isolated vertex in this forest and the remaining six vertices constitute atree. Similarly, the subgraph induced by the vertices colored blue or green is a forest
consisting of the tree on the six vertices a4 , a6 , a7 , a8 , a9 , a10 and the isolated vertex
a2 ; and the subgraph induced by the vertices colored red or green is a forest consistingof the path on the four vertices a3 , a6 , a5 , a8 and the two isolated vertices a1 and a10 .
Our goal here is to show how the ordering in a substitutable partition can becomputed from an acyclic coloring of the adjacency graph G(A) of the symmetricmatrix A. First, it is easy to verify that in a group of columns in a substitutablepartition induced by an acyclic coloring, each diagonal element is the only nonzero
in its row. This is an immediate consequence of the fact that an acyclic coloring is
a distance-1 coloring, for if aij is nonzero, vertices ai and aj are adjacent and areassigned diﬀerent colors in a distance-1 coloring. Thus ai is the only column with anonzero in row i in the group of columns to which it belongs.
Next, consider oﬀ-diagonal elements aij with i = j, where the two elements aij
and aji are identiﬁed due to symmetry. Each edge in G(A) belongs to exactly one twocolored subgraph of G(A), since its two endpoints are assigned distinct colors. Thusan acyclic coloring partitions the edges of the graph G(A) into two-colored subgraphs,
which are forests. We will solve for the matrix elements corresponding to the edges
from each two-colored tree separately, and thus we need to prescribe the order only
for a tree.
In a two-colored tree, every edge incident on a leaf vertex can be solved for
immediately since it is the only nonzero in a row of the group of columns to which
its parent vertex belongs. As an example, consider the red-blue forest induced by
the vertices {a1 , a2 , a3 , a4 , a5 , a7 , a9 } in Figure 6.1. The relevant nonzeros in the twogroups of columns {a1 , a3 , a5 } and {a2 , a4 , a7 , a9 } in the partition induced by theacyclic coloring are as follows:


a21





 a71

a23
a43
a63

a25





a65 
;

a85



a12
 a32

 a52





a17
a34

a87
a10,4





.
a69 

a89 
a10,9

GRAPH COLORING FOR COMPUTING DERIVATIVES

673

Note that we have not shown the diagonal nonzeros in each column. Furthermore,
the elements shown in bold will be evaluated when we consider the other two-coloredforests; the bold entries in the submatrix on the left will be evaluated from the redgreen forest, and the bold entries in the submatrix on the right will be evaluated fromthe blue-green forest.
Consider the nontrivial tree in the red-blue forest. In this tree, the vertex a7 is
a leaf, and the edge (a7 , a1 ) can be computed from column a1 since a71 is the only
nonzero in row 7 in the columns colored red, the set {a1 , a3 , a5 }. The vertex a4 is alsoa leaf, and the edge (a4 , a3 ) can be computed from column a3 since a43 is the only
nonzero in row 4 in the columns colored red. Likewise, vertex a5 is a leaf, and theedge (a5 , a2 ) can be computed from column a2 since a52 is the only nonzero in row 5in the columns colored blue, the set {a2 , a4 , a7 , a9 }.
Furthermore, once the edges incident on a leaf have been evaluated, they canbe deleted from the tree to create new leaves. The process can be repeated to solvefor edges incident on the new leaf vertices, by using values computed for the leafedges from earlier steps. The deletion of edges corresponds to substitution of alreadycomputed quantities in a triangular system of equations. In our illustration, once edges(a7 , a1 ), (a4 , a3 ), and (a5 , a2 ) have been evaluated and deleted, the path a1 , a2 , a3remains. In this path, vertices a1 and a3 are leaves. In the columns colored blue,row 1 has two nonzeros, a12 and a17 ; since the latter is known from the evaluationof a71 , we can compute the element a12 by substitution. Finally, the last edge inthe tree, edge (a2 , a3 ), can be evaluated by using row 1 of the columns colored red,by substituting the known values for a21 and a25 , or by using row 2 of the columns
colored blue, by substituting the known value for a34 .
To summarize, an ordering from which matrix elements can be evaluated follows
naturally from an acyclic coloring by considering the partition of the graph into twocolored trees. In each tree, the edges incident on a leaf vertex can be evaluateddirectly; subsequently, as the evaluated edges are deleted, edges incident on new leaf
vertices could be evaluated in a substitution method. This process can be repeateduntil the tree becomes empty and every edge in the tree has been evaluated.
Coleman and Cai [27] established the connection between acyclic coloring and thecomputation of a symmetric matrix using a substitution method. Acyclic coloring had
been studied earlier by Grünbaum [56] in a diﬀerent context. Further information on
acyclic coloring with pointers to references is included in section 11.4.
Theorem 6.4 (Coleman and Cai [27]). Let A be a symmetric matrix with nonzero
diagonal elements and let G(A) = (V, E) be its adjacency graph representation. A
mapping φ is an acyclic coloring of G(A) if and only if φ induces a substitutablepartition of the columns of A.
Proof. Assume that φ is an acyclic coloring of the graph G(A). We claim thatφ induces a substitutable partition of the columns of A. This is immediate fromthe earlier discussion, since an acyclic coloring partitions the edges of G(A) into twocolored forests, and the edges in each two-colored tree can be computed by the processof computing edges incident on leaf vertices.
Now we prove the converse. Assume that we are given a substitutable partitionof the columns C = {C1 , C2 , . . . , Cp } of the symmetric matrix A. Deﬁne a mapping
φ(ai ) = k if ai ∈ Ck for 1 ≤ i ≤ p. We will show that φ corresponds to an acyclic
coloring of the vertices of the adjacency graph G(A).
To arrive at a contradiction, assume that φ is not an acyclic coloring. Then either
it is not a distance-1 coloring, or the graph G(A) has a cycle in which the vertices are
assigned only two colors.

674

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

If φ is not a distance-1 coloring, then there are two adjacent vertices, say, ai andaj in the graph G(A) with the same color. Then these columns belong to the samegroup in C. The columns ai and aj contain the following 2 × 2 submatrix:
aii aij.
aji ajj
These three distinct nonzero elements are not present in any other column and hencecannot be evaluated from any other group. But their values cannot be computed bysolving a triangular system of equations since there are two unknown elements in eachof the two rows. This contradicts the assumption that C is a substitutable partition.
Now consider the case when φ induces a two-colored cycle in the graph G(A). By
assumption, the elements corresponding to the edges of the cycle can be computedfrom the substitutable partition of A. Let (ai , aj ) denote the ﬁrst edge computed fromthe two-colored cycle. Let ah denote the neighbor of ai distinct from aj in the cycle,
and similarly let ak denote the neighbor of aj distinct from ai in the cycle. Then,by choice of the edge (ai , aj ), the edges (ah , ai ) and (aj , ak ) have not been evaluatedyet. Thus the pair of nonzeros ahi and aij (and also aij and ajk ) cannot be ordered
to satisfy the requirements of a substitutable partition given in Deﬁnition 6.1. This
contradicts the assumption that C is a substitutable partition.
In view of Theorem 6.4, Problem 6.2 is equivalent to the following graph problem.
Problem 6.5. Given the adjacency graph G(A) = (V, E) representing the sparsity structure of a symmetric matrix A with nonzero diagonal elements, ﬁnd an acyclic
coloring of G(A) that uses the fewest colors.
6.1.3. An Acyclic Coloring Algorithm. Coleman and Cai [27] proved that the
problem of ﬁnding an acyclic coloring of a graph that uses the fewest colors is NP-hard.
Recently, we developed a new eﬃcient heuristic algorithm for the problem. Here webrieﬂy mention the main idea in the algorithm; for a more detailed discussion see [48].
Our algorithm relies on the fact that in a graph where the vertices have beencolored satisfying the conditions of acyclic coloring, every subgraph induced by a setof vertices assigned any two colors is a forest. The union of all two-colored forests
in a partially acyclically colored graph is a collection of edge-disjoint trees. Eachtree in the collection is dynamic in the sense that edges may be added to it as thealgorithm proceeds. Based on these observations, our algorithm uses the disjoint-set
data structure to maintain the collection of trees in an eﬃcient way. In our context, a
set in the disjoint-set data structure corresponds to a two-colored tree and an elementof a set corresponds to an edge. The two important operations supported by the
disjoint-set data structure are union and ﬁnd.
The algorithm iterates over the vertices in the set V in some order. In each stepi, 1 ≤ i ≤ |V |, it ensures that the color chosen for vertex vi (1) is distinct from eachof the colors of the distance-1 neighbors of vi , and (2) does not lead to a cycle inany one of the current two-colored trees. The ﬁrst requirement is easy to enforce:the colors used by the distance-1 neighbors of vertex vi are excluded from the setof colors allowed for vertex vi . Since we maintain a collection of two-colored trees,
the second requirement reduces to checking whether vertex vi is connected to at least
two same-colored vertices in a single two-colored tree. These checks can be done
systematically for all the trees incident on vi in order to exclude additional colors
from the set of colors allowed to vertex vi . Once the colors that are not allowed forvertex vi are determined, the smallest allowable color is chosen and assigned to vi . The
computational work involved in coloring the vertex vi is proportional to d2 (vi ). The

GRAPH COLORING FOR COMPUTING DERIVATIVES

675


total number of ﬁnd operations in the algorithm is bounded by f = v∈V d2 (v), while
the number of union operations is at most |E|. Thus the overall time complexity of
the algorithm is O(α(f, |E|) · |V |δ 2 ), where α is the functional inverse of Ackermann’sfunction. Recently, we used the underlying technique in this algorithm, exploiting thestructure of two-colored induced subgraphs, to design an O(|V |δ 2 )-time star coloringalgorithm. See the report [48] for a detailed discussion of these algorithms.
Earlier, Coleman and Cai [27] suggested an algorithm for the acyclic coloring
problem. The idea in their algorithm is to ﬁrst transform a given graph G = (V, E) toa “completed” graph G = (V, E  ) by adding edges in such a way that every cycle in
G includes a triangle in G , and then use a known distance-1 coloring heuristic on G .
The construction of G is done in the following way. Initially set E  to be the same as
E. Consider an ordering (numbering) π : V ↔ {1, 2, . . . , n} of the vertices. Process
the vertices in the order given by π. In each iteration i, if the vertex w = π −1 (i)is such that there exists a path v, w, x in G where π(w) > max{π(v), π(x)}, then
add the edge (v, x) to E  . It is clear that when G is constructed in this manner,a distance-1 coloring of G gives a valid acyclic coloring of G. However, an optimal
distance-1 coloring of G does not imply an optimal acyclic coloring of G. In fact, thealgorithm of Coleman and Cai is an algorithm for a more constrained coloring problem
that Coleman and Moré [31] have called triangular coloring. (A vertex coloring φ ofa graph G is triangular if there exists an ordering π of the vertices such that (1) φis a distance-1 coloring of G, and (2) in every path v, w, x in G, the vertices v and x
receive diﬀerent colors whenever π(w) > max{π(v), π(x)}.)The graph G used in the algorithm described above depends critically on thechosen ordering π. Coleman and Moré [31] showed that a smallest last vertex ordering
is particularly well-suited for minimizing the number of colors required to distance-1
color the graph G . For a discussion of the smallest last ordering and its relationshipto various graph parameters, see section 11.
The graph G has interesting relationships with ﬁll graphs in sparse matrix factorization. Our ﬁrst observation is that the graph G is a subgraph of a chordal
completion G+ = (V, E ∪ F ) of G. A graph is chordal if every cycle on four or morevertices has a chord, an edge connecting two nonconsecutive vertices on the cycle.
Computing a chordal completion of a graph in which the number of ﬁll edges |F | isminimized is NP-hard [124]. This graph problem is a well-known model for sparse
Cholesky factorization. The process through which G is computed resembles, but is
not the same as, the elimination game, a process through which G+ is computed. The
diﬀerence is that in the former case a ﬁll edge is added between two lower numberedvertices whenever these vertices have a common higher numbered neighbor in theoriginal graph G.
We can also characterize the edges in the graph G in terms of ﬁll in a sparse
incomplete LU factorization. In a level-based incomplete LU factorization, each of
the edges in the original graph are assigned level 0, and the ﬁll edges are assigned
higher values of levels. In the sum rule of assigning levels to ﬁll edges, when a ﬁll edge
(vj , vk ) is created by a pair of edges (vi , vj ) and (vi , vk ), the level of the ﬁll edge is
one more than the sum of the levels of the two causative edges. In this case, since thecausative edges belong to the original graph G, they have level values equal to 0, andhence the ﬁll edge (vj , vk ) has level 1. Thus the added edges E  \ E in G are exactly
the ﬁll edges in a level-based sparse incomplete LU factorization corresponding tolevel 1, i.e., ILU(1).
Hysom and Pothen [70] proved a theorem characterizing ﬁll in an incomplete
factorization. A ﬁll path between two vertices vi and vj is a path in the original graph

676

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

Algorithm 6.1. A scheme for solving the acyclic bicoloring problem.
procedure AcyclicBicoloringScheme(Gb = (V1 , V2 , E))
1. Find a suitable vertex cover C in Gb
2. Assign the vertices in the set I = (V1 ∪ V2 ) \ C the color 0
3. Color the vertices in C such that the result is an acyclic bicoloring of Gb
end procedure

G joining the two vertices such that all of its interior vertices are numbered less thanvi and vj . An edge (vi , vj ) is a ﬁll edge of level if and only if the shortest ﬁll pathjoining vi and vj in the graph G has length + 1 edges [70].
6.2. Computing the Jacobian. In estimating a nonsymmetric matrix using asubstitution method, the requirement on the bidirectional partition can be relaxed so
as to obtain fewer groups compared to a partition in a direct method.
6.2.1. The Matrix Partitioning Problem. We state the following deﬁnition, due
to Coleman and Verma [32], to subsequently describe the ﬁfth matrix partitioning
problem of our concern.
Definition 6.6. A bidirectional partition (ΠC , ΠR ) of a matrix A is substitutable if there exists an ordering on the elements of A such that for every nonzero
element aij , either (1) column aj is in a group where all nonzeros in row ri , fromother columns in the group, are ordered before aij , or (2) row ri is in a group where
all the nonzeros in column aj , from other rows in the group, are ordered before aij .
Problem 6.7. Given the sparsity structure of an m × n matrix A, ﬁnd a substitutable bidirectional partition (ΠC , ΠR ) of A such that |ΠC | + |ΠR | is minimized.
6.2.2. A Graph Coloring Formulation. The relationship between bicoloring andbidirectional partition, established by Theorem 5.6, coupled with that between acyclic
coloring and substitutable partition, established by Theorem 6.4, suggests that “acyclic
bicoloring” might be the right graph model for Problem 6.7. Coleman and Verma [32]showed that this was indeed the case.
Definition 6.8. Let Gb = (V1 , V2 , E) be a bipartite graph. A mapping φ :[V1 , V2 ] → {0, 1, . . . , p} is an acyclic bicoloring of Gb if the following conditions hold:
1. If u ∈ V1 and v ∈ V2 , then φ(u) = φ(v) or φ(u) = φ(v) = 0.
2. If (u, v) ∈ E, then φ(u) = 0 or φ(v) = 0.
3. If vertices u and v are adjacent to a vertex w with φ(w) = 0, then φ(u) = φ(v).
4. Every cycle uses at least three colors.
Theorem 6.9 (Coleman and Verma [32]). 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
bicoloring of Gb (A) if and only if φ induces a substitutable bidirectional partition(ΠC , ΠR ) of A.
By Theorem 6.9, the following coloring problem is equivalent to Problem 6.7.
Problem 6.10. Given the bipartite graph Gb (A) = (V1 , V2 , E) representing thesparsity structure of an m × n matrix A, ﬁnd an acyclic bicoloring of Gb (A) that uses
the fewest colors.
6.2.3. An Acyclic Bicoloring Algorithm. Observation 5.8 suggests the scheme
AcyclicBicoloringScheme outlined in Algorithm 6.1 for solving the acyclic bicoloring problem. Step 3 of AcyclicBicoloringScheme can be done using anappropriate adaptation of the acyclic coloring algorithm sketched in section 6.1.3.

GRAPH COLORING FOR COMPUTING DERIVATIVES

677

7. Interrelationships among the Coloring Problems. In this short section, we
expose the interrelationships among the various coloring variants introduced thus far
in this paper. The relationships reveal that distance-2 coloring is the most generalvariant.
7.1. Chromatic Numbers. The conditions required by distance-1 coloring, acyclic coloring, star coloring, acyclic bicoloring, star bicoloring, and distance-2 coloring
imply the relationships among the respective chromatic numbers stated in Theorems 7.1 and 7.2. The two theorems are due to Coleman and Cai [27] and Coleman
and Verma [32].
Recall that the distance-k chromatic number of a graph G is denoted by χk (G),and the least number of colors required for a partial distance-2 coloring of a bipartite
graph Gb = (V1 , V2 , E) on Vi (i = 1, 2) is denoted by χ2 (Gb , Vi ). Let the chromatic
number for acyclic and star coloring of a graph G be denoted by χa (G) and χs (G),respectively. Further, let the chromatic number for acyclic and star bicoloring of abipartite graph Gb be denoted by χab (Gb ) and χsb (Gb ), respectively.
Theorem 7.1. For every graph G = (V, E),
χ1 (G) ≤ χa (G) ≤ χs (G) ≤ χ2 (G) = χ1 (G2 ).
Proof. Regarding the ﬁrst three inequalities, observe that a distance-2 coloringis a star coloring; a star coloring is an acyclic coloring; and an acyclic coloring is
a distance-1 coloring. The last equality holds since a distance-2 coloring of G isequivalent to a distance-1 coloring of G2 .
Theorem 7.2. For every bipartite graph Gb = (V1 , V2 , E),
χab (Gb ) ≤ χsb (Gb ) ≤ min{χ2 (Gb , V1 ), χ2 (Gb , V2 )}.
Proof. The ﬁrst inequality is obvious. For the second inequality, observe that a
partial distance-2 coloring on V2 is a valid star bicoloring of the bipartite graph Gb
where all the vertices in V1 are restricted to be colored with color zero. A similarargument, with the roles of V1 and V2 interchanged, can be used to complete theproof.
In the context of eﬃcient derivative matrix computation, the implication of Theorem 7.2 is that an optimal bidirectional partition, irrespective of the structure of thematrix, yields at most as many groups as an optimal unidirectional partition, andhence potentially results in a more eﬃcient computation.
Theorems 7.1 and 7.2 show that distance-2 coloring is an archetypal model in thecomputation of Jacobian and Hessian matrices.
Distance-2 coloring also has applications other than derivative matrix computation. Examples include channel assignment [84] and facility location problems (seeChapter 5 of the book by Vazirani [117]).
7.2. Two-Colored Induced Subgraphs. As discussed in earlier sections, thenames for star and acyclic coloring are derived from the structure of two-coloredinduced subgraphs in the respective colorings. In this section we revisit this issue andextend the observation to the distance-1 and distance-2 coloring cases.
Let G = (V, E) be a graph and φ : V → {1, 2, . . . , p} be a variant of vertex
coloring. Let H2colors ⊆ G denote a subgraph of G induced by the union of any
two color classes in φ. In Table 7.1 we summarize our observations regarding thecharacterization of H2colors in the cases where φ is a distance-1, an acyclic, a star,
and a distance-2 coloring. The reader will ﬁnd the illustration in Figure 7.1 helpful

678

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN
Table 7.1 Characterization of the subgraph induced by the union of any two color classes.

φDistance-1 coloring
Acyclic coloring
Star coloring
Distance-2 coloring
d

b

d

g

b

e

d

g

b

e

a

i

c

H2colors
A bipartite graph
A collection of trees
A collection of stars
A distance-2 matching

i

c

g

b

e

a

h

d

h

e

a

i

c

g

h

a

i

c

h

f

f

f

f

i

ii

iii

iv

Fig. 7.1 (Left to right) A distance-1, an acyclic, a star, and a distance-2 coloring of a graph.
Table 7.2 The color classes in each of the four colorings depicted in Figure 7.1.

Class 1
Class 2
Class 3
Class 4
Class 5Class 6i{a, d, e, f, i}{b, c}{g, h}

ii{a, d, f, i}{b, c}{e}{g, h}

iii{a, i}{b, c}{d, f }
{e}{g, h}

iv{a, i}{b, f }{c, d}{e}{g}{h}in following this discussion; parts (i) through (iv) depict an example of each of therespective coloring variants. For convenience, Table 7.2 lists the color classes in eachof the four cases of Figure 7.1.
It is obvious that φ is a distance-1 coloring if and only if H2colors is a (notnecessarily connected) bipartite graph.
Since φ in Table 7.1 gets progressively restricted as one goes down the column, correspondingly, H2colors becomes progressively restricted (and less connected). Acyclic
coloring is named as such precisely because the corresponding H2colors is a collectionof trees (an acyclic graph). Similarly, star coloring owes its name to the fact that thecorresponding H2colors is a collection of stars [41]. It is easy to see that a connected
two-colored induced subgraph here ought to be a star, for otherwise there would exist
a two-colored path on four vertices, violating a condition of star coloring.
In the case where φ is a distance-2 coloring, we note that H2colors is a restricted
matching. In the usual sense, a matching in a given graph G = (V, E) is a set M ⊆ E
of edges with no shared endpoints. In other words, in a connected graph G, a pairof edges in a matching M is at least one edge apart in G. The matching H2colors in
the last row of Table 7.1 is such that a pair of edges in H2colors is apart by a pathof length at least two edges in G. We call this a distance-2 matching. The subgraphinduced by the vertices {b, c, d, f } in part (iv) is an example of a distance-2 matching,since it consists of the two edges (b, d) and (c, f ), which are apart from each other bypaths of length two or greater.

GRAPH COLORING FOR COMPUTING DERIVATIVES

679

As a corollary, we observe a relationship between distance-2 vertex coloring anda restricted variant of edge coloring. In the standard usage, an edge coloring of agraph is an assignment of colors to its edges such that every pair of edges sharing an
endpoint receives diﬀerent colors. Such an edge coloring clearly partitions the edge
set into matchings. A distance-2 edge coloring is an assignment of colors to edges suchthat every pair of edges that either shares an endpoint or is separated by one otheredge is assigned diﬀerent colors. Consider now the edge coloring φ derived from a
distance-2 vertex coloring φ in the following way: assign each edge a color obtainedby “mixing” the colors of its two endpoints; i.e., for each e = (u, v), φ (e) = φ(u) 
φ(v). Here, x  y denotes the number obtained by concatenating the positive integers
x and y; e.g., 2  3 is 23. It is evident that the edge coloring φ obtained in thismanner is a distance-2 edge coloring and hence partitions the edge set into distance-2
matchings.
8. Partial Matrix Computation. In many large-scale optimization contexts, theJacobian or the Hessian is formed only for preconditioning purposes, and only a subsetof the matrix elements needs to be computed. Computing a good preconditioner iscritical for fast convergence to a solution. A recent survey article by Knoll and Keyes[79] discusses various applications where the “Jacobian-free Newton–Krylov” method
is used. A basic ingredient of this method is an approximate computation of someelements of the Jacobian. Also, there are examples in which only certain elements ofthe Hessian need to be updated in an iterative procedure, since the other elements do
not change in value [11].
In this section we develop graph coloring formulations of partitioning problems
that arise in the computation of a speciﬁed subset of the nonzero entries of a matrix.
We call this partial matrix computation as opposed to full matrix computation, thecase in which all nonzero entries are required to be determined.
The coloring formulations in this section are new and more sophisticated than thecoloring formulations in full matrix computation. In a coloring formulation of partial
matrix estimation, one may be tempted to think that the vertices corresponding to thenonzeros outside the set of required elements might as well be ignored. However, suchvertices still need to be considered since they could interfere with the computation of
the required elements. The motivation for developing the new graph formulations isthat eﬃcient partial matrix computation can be used to further reduce the numberof colors needed to compute the required elements. For example, if only the diagonal
elements of a Hessian are needed, then we need only a distance-1 coloring of theadjacency graph, rather than a star coloring, the coloring required for full matrix
computation.
Similar to the bicoloring cases discussed in earlier sections, in the colorings deﬁned
here we allow a vertex to be assigned the “neutral” color zero. A vertex with color
zero signiﬁes the fact that a column or a row that corresponds to the vertex is not
used to compute any element in that column or row.
The rest of this section is organized in three parts. Each part deals with a scenariodeﬁned by the kind of matrix under consideration (Jacobian or Hessian) and the typeof partition employed (unidirectional or bidirectional). In each case, the required
entries are assumed to be determined using a direct method. The problems that
correspond to computation via substitution are not studied in this work.
8.1. Unidirectional Computation of the Jacobian. Let A be an m × n nonsymmetric matrix, and let S denote the set of nonzero elements of A required to becomputed. A partition of a subset of the columns of A is structurally orthogonal when680

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

1 2 3 4 5 6
X
1 X
X X
2 X
X
X
3
X
X X
4
X X
5
X X
6

Fig. 8.1

r1
r2
r3
r4
r5
r6

a1
a2
a3
a4
a5
a6

Partial estimation of a nonsymmetric matrix and its formulation as a restricted distance-2
coloring in the associated bipartite graph. The required matrix entries are encircled andthe corresponding edges in the graph are shown in bold. The illustration uses two colors:
columns a1 and a5 are green, and columns a4 and a6 are red.

restricted to S if for every aij ∈ S, column aj is included in some group that contains
no other column with a nonzero in row ri . (The latter nonzero could be inside oroutside S.) Such a partition enables a direct determination of the elements of S.
Let Gb (A) = (V1 , V2 , E) be the bipartite graph of A, and let ES ⊆ E correspondto the elements of S. A mapping φ : V2 → {0, 1, . . . , p} is a distance-2 coloring of
Gb when restricted to ES if the following conditions hold for every edge (v, w) ∈ ES ,where v ∈ V1 , w ∈ V2 :1. φ(w) = 0, and
2. for every path (u, v, w), φ(u) = φ(w). (Here the edge (u, v) may or may notbelong to ES .)Theorem 8.1 below shows the equivalence between the restricted partitioning andcoloring variants introduced above. See Figure 8.1 for an illustration; in the ﬁgure, thenonzero matrix entries required to be computed are encircled and the correspondingedges in the bipartite graph are shown in bold; columns a1 and a5 are colored green,and columns a4 and a6 are colored red.
Theorem 8.1. The mapping φ is a distance-2 coloring of Gb (A) when restricted
to ES if and only if φ induces a structurally orthogonal column partition when restricted to S.
Proof. Assume that φ is a distance-2 coloring of Gb (A) when restricted to ES .
Let p be the number of colors used. We show that the groups {C1 , . . . , Cp } whereCα = {aj : φ(aj ) = α}, 1 ≤ α ≤ p, constitute a structurally orthogonal column
partition when restricted to S. First, observe that by condition 1, for every aij ∈ S(i.e., (ri , aj ) ∈ ES ), φ(aj ) = 0. Thus column aj belongs to group Cφ(aj ) and hence
is involved in the partition. Assume now that the partition induced by the coloring
is not structurally orthogonal when restricted to S. This occurs only if there exist
nonzero elements aij and aik , j = k, such that aij ∈ S and both aj and ak belong
to group Cα for some α , 1 ≤ α ≤ p. But this contradicts condition 2, and hencecannot occur.
Conversely, assume that the partition C = {C1 , . . . , Cp } is structurally orthogonal
when restricted to S. Construct a coloring φ of Gb (A) as follows: φ(aj ) = α if aj ∈ Cα ,and φ(aj ) = 0 if aj does not belong to any group in C. We claim that φ is a distance2 coloring of Gb (A) when restricted to ES . Each vertex in V2 incident on an edgein ES corresponds to a column with an entry in S, and hence gets a nonzero color.
Thus φ satisﬁes condition 1. Consider any path (aj , ri , ak ) where (ri , aj ) ∈ ES . Notethat such a path in Gb (A) exists whenever entries aij and aik are nonzero. Structural

681

GRAPH COLORING FOR COMPUTING DERIVATIVES

1 2 3 4 5 6
X X

X X

X X X

X X

a1

X X X X
X

X X
X X X
X X

Fig. 8.2

X
X

X

a5

a2
a6

a3
a4

Partial estimation of a symmetric matrix and its formulation as a restricted star coloring in the adjacency graph. The required matrix entries are encircled, and edges of thegraph corresponding to required oﬀ-diagonal matrix entries and vertices “corresponding”
to required diagonal matrix entries are shown in bold. The illustration uses three colors:
column a2 is red, a3 is green, and both a4 and a5 are blue.

orthogonality when restricted to S implies that column ak cannot be in the same groupas column aj . Thus, by construction, φ(aj ) = φ(ak ), satisfying condition 2.
8.2. Computing the Hessian. Let A be a symmetric matrix with nonzero diagonal elements, and let S denote the set of nonzero elements of A required to becomputed. A partition of a subset of the columns of A is symmetrically orthogonal
when restricted to S if for every aij ∈ S at least one of the following two partition
conditions are met:
1. The group containing column aj has no other column with a nonzero in row
ri .
2. The group containing column ai has no other column with a nonzero in row
rj .
Such a partition enables a direct determination of the elements of S.
Let G(A) = (V, E) be the adjacency graph of A; let Sod ⊆ E correspond tothe oﬀ-diagonal elements in S; and let Sd correspond to the diagonal elements in
S, i.e., Sd = {(u, u) : u ∈ U } where U ⊆ V . Let ES = Sod ∪ Sd . A mapping
φ : V → {0, 1, 2, . . . , p} is a star coloring of G when restricted to ES if the following
coloring conditions hold:1. For every (u, u) ∈ Sd ,1.1. φ(u) = 0, and1.2. for every (u, v) ∈ E, φ(u) = φ(v).
2. For every (v, w) ∈ Sod ,2.1. φ(v) = φ(w), and2.2. at least one of the following two conditions holds:
2.2.1. φ(v) = 0 and for every path (v, w, x), φ(v) = φ(x), or
2.2.2. φ(w) = 0 and for every path (u, v, w), φ(u) = φ(w).
The following theorem states the equivalence between the restricted partitioning
and coloring variants introduced above. See Figure 8.2 for an illustration; in the ﬁgure,the matrix entries required to be computed are encircled, edges in the adjacency graphcorresponding to the required oﬀ-diagonal entries are shown in bold, and column
vertices containing a required diagonal element are shown in bold. In the illustration,column a2 is colored red, a3 is green, and a4 and a5 are blue.
Theorem 8.2. The mapping φ is a star coloring of G(A) when restricted to ES
if and only if φ induces a symmetrically orthogonal partition when restricted to S.

682

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

Proof. Assume that φ is a star coloring of G(A) when restricted to ES . Let
the number of colors used by φ be p. We show that the groups {C1 , . . . , Cp } whereCα = {aj : φ(aj ) = α}, 1 ≤ α ≤ p, constitute a symmetrically orthogonal partition
when restricted to S.
By coloring conditions 1.1 and 2.1, for every aij ∈ S (i.e., (ai , aj ) ∈ ES ), at least
one of the vertices ai or aj has a nonzero color and hence the column is involved in
the partition {C1 , . . . , Cp }. Let aij ∈ S be a diagonal entry (i = j). Then coloring
conditions 1.1 and 1.2 ensure that φ(ai ) = 0 and that φ(ai ) = φ(ak ) for every (ai , ak ) ∈E. Thus, by construction, column ai belongs to group Cφ(ai ) and no column ak , wherei = k and aik = 0, is in Cφ(ai ) . This clearly satisﬁes the partition conditions.
Let aij ∈ S now be an oﬀ-diagonal entry (i = j). Assume without loss of generality that φ(aj ) = 0. By coloring condition 2.1, φ(ai ) = φ(aj ). By coloring condition
2.2.2, there is no path (ah , ai , aj ) in G(A), for any h = i, j, such that φ(ah ) = φ(aj ).
The last two statements together imply that column aj belongs to group Cφ(aj ) andthat no column ah , where h = j and aih = 0, is in Cφ(aj ) . This satisﬁes the ﬁrst
partition condition. A similar argument applies to the case where φ(ai ) = 0, which
implies the satisfaction of the alternate partition condition.
To prove the converse, assume that the partition C = {C1 , . . . , Cp } is symmetrically orthogonal when restricted to S. Construct a coloring φ of G(A) as follows.
Deﬁne φ(aj ) = α if aj ∈ Cα , and φ(aj ) = 0 if aj does not belong to any group in C.
We claim that φ is a star coloring of G(A) when restricted to ES .
Consider a diagonal element aii ∈ S. The partition conditions ensure that column
ai is in some group Cα and that there is no column ak ∈ Cα , i = k, such that aik = 0.
Thus, by construction, φ(ai ) = 0 and φ(ai ) = φ(ak ) for every (ai , ak ) ∈ E, satisfying
coloring condition 1.
Consider now the case where aij ∈ S is an oﬀ-diagonal element. First, observethat since all diagonal elements are nonzero, ai and aj cannot belong to the samegroup. Thus φ(ai ) = φ(aj ), satisfying coloring condition 2.1. Second, observe that
there are two possibilities by which the partitioning conditions can be satisﬁed. We
consider only one of these; the second can be treated in a similar manner. Suppose
column aj belongs to some group Cα and that there is no other column ah ∈ Cα ,h = j, such that aih = 0. Thus, by construction, φ(aj ) = 0 and φ(ah ) = φ(aj ) for
every path (ah , ai , aj ) in G(A), satisfying coloring condition 2.2.2.
A special case of Theorem 8.2 is the problem of estimating only the diagonal
elements of A, i.e., Sd = {(v, v) : v ∈ V } and Sod = ∅. For this problem, coloring
conditions 1.1 and 1.2 are the only applicable conditions, and they imply that a
distance-1 coloring of G(A) suﬃces.
8.3. Bidirectional Computation of the Jacobian. Let A be an m × n nonsymmetric matrix, and let S denote the set of nonzero elements of A required to becomputed. A bidirectional partition (ΠC , ΠR ) of a subset of the columns and rows of
A is structurally orthogonal when restricted to S if for every aij ∈ S at least one of
the following two partition conditions are met:
1. The group containing column aj has no other column with a nonzero in row
ri .
2. The group containing row ri has no other row with a nonzero in column aj .
Such a partition enables a direct determination of the elements of S.
Let Gb (A) = (V1 , V2 , E) be the bipartite graph of A, and let ES ⊆ E correspondto the elements in S. A mapping φ : [V1 , V2 ] → {0, 1, . . . , p} is a star bicoloring of Gb
when restricted to ES if the following coloring conditions are met:

683

GRAPH COLORING FOR COMPUTING DERIVATIVES

1 2 3 4 5
1
2
3
4
5

Fig. 8.3

X

X

X X
X

X
X

X X X X
X

X

r1
r2
r3
r4
r5

a1
a2
a3
a4
a5

Partial estimation of a nonsymmetric matrix via bidirectional partition and its formulation as a restricted star bicoloring in the bipartite graph. The required matrix entries areencircled and the corresponding edges in the graph are shown in bold. The illustration uses
three colors: both r2 and r3 are blue, r4 is green, and a1 is red.

1. Vertices in V1 and V2 receive disjoint colors, except for color 0;
i.e., for every u ∈ V1 and v ∈ V2 , either φ(u) = φ(v) or φ(u) = φ(v) = 0.
2. At least one endpoint of an edge in ES receives a nonzero color;i.e., for every (v, w) ∈ ES , φ(v) = 0 or φ(w) = 0.
3. For every edge (v, w) ∈ ES ,3.1. if φ(v) = 0, then, for every path (u, v, w), φ(u) = φ(w);3.2. if φ(w) = 0, then, for every path (v, w, x), φ(v) = φ(x);3.3. if φ(v) = 0 and φ(w) = 0, then for every path (u, v, w, x), either φ(u) =φ(w) or φ(v) = φ(x).
The following theorem establishes the equivalence between the restricted partitioning and coloring notions introduced above. See Figure 8.3 for an illustration.
Theorem 8.3. The mapping φ is a star bicoloring of Gb when restricted to ES if
and only if φ induces a structurally orthogonal bidirectional partition when restricted
to S.
Proof. Let the construction of a partition given a coloring, and vice versa, bedone in a similar manner as in the proof of Theorem 8.1.
Assume that φ is a star bicoloring of Gb (A) when restricted to ES . Let the induced
bidirectional partition be (ΠC , ΠR ). Coloring condition 1 implies that (ΠC , ΠR ) is abidirectional partition. By condition 2, for every aij ∈ S, either aj ∈ ΠC or ri ∈ ΠR(or both). Assume now that (ΠC , ΠR ) is not structurally orthogonal when restricted
to S. This occurs only if one of the following cases holds for any nonzero aij ∈ S:
• φ(ri ) = 0, φ(aj ) = 0, and there exists a column ak , where j = k and aik = 0,such that φ(aj ) = φ(ak ). But this contradicts coloring condition 3.1 andhence cannot occur.
• φ(aj ) = 0, φ(ri ) = 0, and there exists a row rh , where h = i and ahj = 0,such that φ(rh ) = φ(ri ). But this contradicts coloring condition 3.2 and hencecannot occur.
• φ(ri ) = 0, φ(aj ) = 0, and there exists a column ak , where j = k and aik = 0,
and a row rh , where h = i and ahj = 0, such that φ(aj ) = φ(ak ) andφ(rh ) = φ(ri ). But this contradicts coloring condition 3.3 and hence cannot
occur.
Hence, the bipartition (ΠC , ΠR ) is structurally orthogonal when restricted to S.
Conversely, assume that (ΠC , ΠR ) is a structurally orthogonal bidirectional partition when restricted to S. Clearly, the constructed coloring φ satisﬁes conditions 1

684

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

and 2. To complete the proof, we show that φ also satisﬁes condition 3. Assume thatφ violates condition 3. Then one of the following cases must hold:• There exists a path (ah , ri , aj ) for some (ri , aj ) ∈ ES such that φ(ri ) = 0 andφ(ah ) = φ(aj ). But this implies that element aij cannot be determined directly, contradicting the assumption that (ΠC , ΠR ) is structurally orthogonal
when restricted to S.
• There exists a path (ri , aj , rk ) for some (ri , aj ) ∈ ES such that φ(aj ) = 0
and φ(ri ) = φ(rk ). Again this implies that element aij cannot be determined
directly, a contradiction of our assumption.
• There exists a path (ah , ri , aj , rk ) for some (ri , aj ) ∈ ES such that φ(ri ) =φ(rk ) = 0 and φ(ah ) = φ(aj ) = 0. But this implies that element aij cannot
be determined directly, contradicting the assumption.
Algorithms. We have not developed specialized algorithms for the restricted coloring
problems introduced in this section. However, we believe that the ideas used in thealgorithms for the corresponding coloring problems in full matrix computation can beadapted to the restricted cases.
9. Hypergraph Coloring Formulations. In this section, we introduce yet another perspective concerning Problem 3.3, the partitioning problem that arises in theunidirectional computation of a Jacobian via a direct method. This view uses thenotion of a hypergraph, which is a generalization of a graph. Besides being an interesting alternative perspective, a hypergraph formulation might be a worthwhile
approach for modeling partial matrix computation problems for preconditioning purposes. However, hypergraph formulations are included here mainly for the sake ofcompleteness.
9.1. Deﬁnitions. We begin by deﬁning a few concepts that we need for our formulations. A hypergraph H = (V, E) consists of a ﬁnite set V of vertices and acollection E of nonempty subsets of V called hyperedges. Note that the size of a hyperedge (i.e., the number of vertices in it) is not constrained to be two as in a graph.
A hypergraph is called r-uniform if all of its hyperedges are of size r. Thus, a graphis a 2-uniform hypergraph.
We denote the vertex set and edge set of a hypergraph H by V (H) and E(H),respectively. When the hypergraph under consideration is clear from the context, we
may use just V and E to refer to these sets.
The line graph L(H) of hypergraph H is a graph where the vertices of L(H) are
the hyperedges of H and the edges of L(H) are pairs of intersecting hyperedges of H.
Figure 9.1 (left) is a simple example that shows a hypergraph consisting of ﬁve
vertices and three hyperedges; two hyperedges have size three and one hyperedge hassize two. The ﬁgure on the right shows the corresponding line graph.
9.2. Hypergraph Coloring. Given a hypergraph H = (V, E), a strong vertex
coloring of H is a mapping φ : V → {1, 2, . . . , p} such that for every hyperedgee ∈ E and every pair of vertices {u, v} ⊆ e, φ(u) = φ(v). In words, we requirethat the vertices in every hyperedge be assigned distinct colors in a strong vertex
coloring of a hypergraph. The more common variant of hypergraph coloring, called
weak vertex coloring, requires that no hyperedge be monochromatic; i.e., at least two
diﬀerent colors are assigned to the vertices in every hyperedge consisting of two or
more vertices.
An edge coloring of a hypergraph H = (V, E) is a mapping φ : E → {1, 2, . . . , p}such that hyperedges e1 and e2 in E satisfy φ(e1 ) = φ(e2 ) whenever e1 ∩ e2 = ∅.

GRAPH COLORING FOR COMPUTING DERIVATIVES

e1
v1

v2

v5

e2

e1

685

e2

v3
v4
e3

e3

Fig. 9.1 A hypergraph and its line graph.

The strong (weak) vertex coloring problem of a hypergraph asks for a strong(weak) vertex coloring with the fewest colors. Similarly, the edge coloring problem ofa hypergraph asks for an edge coloring with the fewest colors.
For more information on the concepts introduced here and other notions in hypergraphs as well as pointers to relevant literature, see, e.g., the book by Berge [13] andthe survey article by Duchet [38]. The terms strong and weak coloring in hypergraphs
are used, for instance, in [83, 112].
9.3. Representing Matrices Using Hypergraphs. To achieve our goal of formulating Problem 3.3 using hypergraphs, we introduce two new ways in which thesparsity structure of a matrix can be represented.
Let A be an m × n matrix. Deﬁne the column-oriented hypergraph representationof A as the hypergraph Hc (A), where V (Hc (A)) = {v1 , v2 , . . . , vn } and each vj corresponds to column aj ; and E(Hc (A)) = {e1 , e2 , . . . , em } with each hyperedge ei = {vj :aij = 0} “corresponding” to row ri . The quotation marks are used to reﬂect that a
row is perceived as a set of columns at which it has nonzero entries. Similarly, deﬁne
the row-oriented hypergraph representation of A as the hypergraph Hr (A), where} and each vi corresponds to row ri ; and E(Hr (A)) =V (Hr (A)) = {v1 , v2 , . . . , vm{e1 , e2 , . . . , en } with each ej = {vi : aij = 0} “corresponding” to column aj .
9.4. Structurally Orthogonal Partition and Hypergraph Coloring. Our main
result in this section is Theorem 9.1, which summarizes the equivalence relationshipsthat underlie the graph- and hypergraph-theoretic formulations of Problem 3.3.
Theorem 9.1. Let A be an m × n matrix. Let Hc and Hr be the column-oriented
and row-oriented hypergraph representations of A, respectively. Let Gb and Gc be thebipartite and column intersection graph representations of A, respectively. Let φ be amapping S → {1, 2, . . . , p}, where S is a context-dependent set of n elements and p isa positive integer satisfying p ≤ n. Then the following statements are equivalent:
1. φ induces a structurally orthogonal column partition of A.
2. φ is a partial distance-2 coloring of Gb on the column vertex set.
3. φ is a distance-1 coloring of Gc .
4. φ is a strong vertex coloring of Hc .
5. φ is an edge coloring of Hr .
6. φ is a distance-1 coloring of the line graph L(Hr ) of the hypergraph Hr .
Proof. The equivalences among statements 1, 2, and 3 have already been established in section 3. Here we show the equivalences among the remaining statements;
the reader will ﬁnd Figure 9.2 helpful in following the discussion.

686

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

a 11 a 12 0
0 a 15a 21 a 22 000a 31 0
a
0
34 0
a
a
0
0
43 44 a 45
A
a1

a2

a5

Hc
Fig. 9.2

Gb

r1

a3
a4

a1
a2
a3
a4
a5

r1
r2
r3
r4

r4

a1

a2

r2
r3

Hr

a3

a5
a4
Gc

Structurally orthogonal partition and (hyper)graph colorings. In this illustration three colors are used: columns a1 and a3 are green, a2 and a4 are red, and a5 is blue.

First note the relationship between hypergraph Hc and column intersection graph
Gc of A: the columns with a nonzero in row ri form a clique in Gc , whereas they
form a hyperedge in Hc . A similar observation can be made regarding Hr and therow intersection graph of A. With this observation it is obvious that a strong vertex
coloring of Hc is equivalent to a distance-1 coloring of Gc . Next observe that theline graph L(Hr ) of Hr is isomorphic to the column intersection graph Gc . Thus
statements 1, 2, 3, 4, and 6 are equivalent. Further, since the distance-1 neighborsof a vertex v from the set of column vertices of Gb constitute a hyperedge in Hr , itfollows that statements 2 and 5 are equivalent, thus completing our proof.
As we pointed out earlier in this section, the hypergraph formulation has a potential application in partial matrix computation for preconditioning. For such purposes,
the idea is to relax the requirement on strong vertex coloring by allowing the reuse ofcolors within a hyperedge. Using such a weaker coloring, one may be able to computea subset of the nonzero entries of a matrix and use the result as a preconditioner.
10. Other Matrix Estimation Methods. The methods for computing Jacobiansand Hessians considered thus far in this paper rely on a number of underlying requirements.
First, each method is based on a unidirectional or a bidirectional partition; i.e.,a set S of columns or rows is divided into disjoint subsets whose union is S. Thisrequirement precludes the possibility of a column (or a row) in S belonging to morethan one group, or to none of the groups. Second, a column or a row of a matrix is seen
as an atomic entity; i.e., a column or a row is not divided into any smaller parts. Third,a partition is required to be structurally orthogonal, or symmetrically orthogonal, or
substitutable. The last requirement implies that the system of equations deﬁned
by the partitions is either diagonal or triangular, and hence the unknowns can beobtained either directly or via substitution.
However, approaches where one or more of these underlying requirements arerelaxed have also been suggested in the literature. In this section, we will discuss

GRAPH COLORING FOR COMPUTING DERIVATIVES

687

a few examples of such approaches. For Hessian estimation, McCormick [97] gives
a classiﬁcation of direct methods, including those that do not necessarily rely on
symmetrically orthogonal partitions.
10.1. Methods Based on Solving a Rectangular System. Newsam and Ramsdell[103] proposed a method for computing a Jacobian that relies on solving an overdetermined rectangular system of equations. The approach enables the determinationof the nonzero entries of an m × n Jacobian matrix using ρmax groups, where ρmax
is the maximum number of nonzeros in a row of the matrix. The method is optimalin terms of the number of groups used. However, it needs to solve n least-squares
problems, and hence the system is potentially ill-conditioned.
Geitner, Utke, and Griewank [49] applied the Newsam and Ramsdell approachwithin the context of automatic diﬀerentiation.
Recently, Hossain and Steihaug [64] suggested an elimination scheme for computing a Jacobian. The scheme is based on successive merging of the columns of thecompressed Jacobian, obtained from a structurally orthogonal partition. The nonzero
entries are then determined by solving a banded system of equations. The latter requires an LU factorization. This approach also uses the optimal number of groups,
the maximum number of nonzeros in a row of the Jacobian matrix.
10.2. Element and Variable Isolation. Newsam and Ramsdell [103] also suggested a method that generalizes structurally orthogonal column partitioning. They
used the term variable isolation (VI) to refer to a method based on a structurally
orthogonal column partition and element isolation (EI) to refer to their generalized
method.
Definition 10.1. A nonzero element aij of a matrix A is isolated from nonzero
element apq , whenever aiq = apj = 0, or j = q and i = p.
The EI method partitions the nonzero entries of a matrix into groups consistingof pairwise isolated elements, with the objective of having the fewest groups. The
nonzeros in a row of a matrix are, by deﬁnition, not isolated from each other; hence,
the maximum number of nonzeros in a row is a lower bound for an optimal EI partition. An EI partition is used to determine a column grouping (in which a group may
contain structurally nonorthogonal columns, and/or columns may belong to several
groups) that enables the determination of all nonzero matrix entries via matrix-vector
products. In particular, given a partition I1 , I2 , . . . , Ip of the nonzero entries of amatrix A into groups of isolated elements, a column grouping C1 , C2 , . . . , Cp is constructed as follows: Ck = {aj : aij ∈ Ik for some i}. In words, a column group Ck
consists of the column indices of the nonzeros in a set Ik .
The middle column in Figure 10.1 illustrates two diﬀerent partitions into isolated elements of a matrix whose nonzero structure is shown on the left. The upper
part shows a partition into the four groups {a11 , a44 , a56 }, {a21 , a34 }, {a13 , a23 , a55 },and {a26 , a35 , a42 }. The column grouping based on this element partition is then{a1 , a4 , a6 }, {a1 , a4 }, {a3 , a5 }, and {a2 , a5 , a6 }. The lower part shows an optimal
partition consisting of three groups: {a11 , a21 , a35 , a42 , a55 }, {a13 , a23 , a34 , a44 }, and{a26 , a56 }. The column grouping corresponding to this EI partition is {a1 , a2 , a5 },{a3 , a4 }, and {a6 }, which in this case coincides with a structurally orthogonal column
partition.
Note that variable isolation in general implies element isolation but not vice versa.
Thus an optimal EI method would use at most as many groups as an optimal VI
method.
Newsam and Ramsdell formulated element isolation as a vertex coloring in agraph. The vertices in the graph are the nonzero entries of the matrix and an edge

688

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

a11 0 a13 0 0 0
a21 0 a23 0 0 a26
0 0 0 a34 a35 00 a42 0 a44 0 0
0 0 0 0 a55 a56

a11 0 a13 0 0 0
a21 0 a23 0 0 a26
0 0 0 a34 a35 00 a42 0 a44 0 0
0 0 0 0 a55 a56

r1
r2
r3
r4
r5

a1
a2
a3
a4
a5
a6

a11 0 a13 0 0 0
a21 0 a23 0 0 a26
0 0 0 a34 a35 00 a42 0 a44 0 0
0 0 0 0 a55 a56

r1
r2
r3
r4
r5

a1
a2
a3
a4
a5
a6

Fig. 10.1 Two partitions of the nonzero entries of a matrix into groups of isolated elements. Eachpartition is also represented as a specialized edge coloring in the bipartite graph.

between two entries exists whenever the latter are not isolated from each other. Sucha coloring formulation is likely to be impractical for large problem instances as it deals
with coloring a graph consisting of nnz vertices and O((nnz)2 ) edges, where nnz isthe number of nonzero entries in the matrix.
10.3. Element Isolation and Edge Coloring. We note that the problem of partitioning the nonzero entries of a matrix into groups of isolated elements can in fact
be modeled as a specialized edge coloring problem on a smaller graph.
Definition 10.2. Let Gb = (V1 , V2 , E) be the bipartite graph of matrix A, where
V1 is the row vertex set and V2 is the column vertex set. A mapping φ : E →{1, 2, . . . , p} is called an EI edge coloring if the following two conditions are satisﬁed:
1. At each row vertex v, edges having v as an endpoint are colored diﬀerently.
2. In every path (e1 , e2 , e3 ) on three edges, φ(e1 ) = φ(e3 ).
From the deﬁnition of isolated elements, it is easy to see that an EI edge coloringin the bipartite graph Gb (A) is equivalent to a partition of the nonzero entries of A
into groups of isolated elements. The right column in Figure 10.1 shows the EI edge
coloring representations of the partitions into isolated elements shown in the middle
column.
10.4. Partitioning Segmented Columns. Hossain and Steihaug [61, 65] suggested a framework that has an EI and a VI method as its special cases. In particular,
they suggest a technique for direct estimation of an m × n Jacobian in which the rows
are ﬁrst grouped into blocks that deﬁne “segmented” columns, and then the segments are partitioned into groups each of which consists of structurally independent
segments. Two segments are said to be structurally independent in a manner analogous to two elements being isolated. The case = 1 corresponds to a VI method andthe case = m corresponds to an EI method. Hossain and Steihaug showed that, forsome matrix structures, there exists an = 1 that results in fewer groups comparedwith an approach based on = 1. In a recent report [65], they characterized thenature of a partition that leads to an optimal number of groups. The optimal partition is also formulated as a vertex coloring in an associated graph. One of the main
results in their work is that any direct determination of a Jacobian via matrix-vector
products implies that the nonzero elements in each product are isolated elements.

689

GRAPH COLORING FOR COMPUTING DERIVATIVES

r1
r2
r3
r4
r5
r6
r7

a11 0 0 a14 0 0
0 a22 0 0 a25 0
0 0 a33 0 0 a36
a41 a42 a43 0 0 0
a51 0 0 0 a55 a56
0 a62 0 a64 0 a66
0 0 a73 a74 a75 0

a1
a2
a3
a4
a5
a6

Fig. 10.2 Eisenstat’s example of a matrix where the nonzeros can be directly determined using fewergroups than are obtained from an optimal structurally orthogonal column partition.

10.5. Eisenstat’s Example. Eisenstat’s (counter)example is cited often in theliterature as a case in point where a Jacobian can be directly determined using fewergroups than that obtained by an optimal structurally orthogonal partition.
The example in its general form is stated for an (n + 1) × n matrix A, wheren
1
is an even integer satisfying n ≥ 6, in the following manner. Matrix A = A
where
A2




n
n
A1 = D1 D2 , A2 = DC3 B0 , and D1 , D2 , D3 ∈ R 2 × 2 are nonsingular diagonalnnmatrices; every diagonal entry of B ∈ R 2 × 2 is zero, while every oﬀ-diagonal entryn
1× n
is nonzero; C ∈ R 2 consists entirely of nonzeros; and 0 ∈ R1× 2 is a zero vector.
Eisenstat’s example for the case n = 6 is shown in Figure 10.2.
The columns in Eisenstat’s example are pairwise structurally nonorthogonal, andhence an optimal structurally orthogonal column partition would require n groups.
However, as has been noted by several authors (e.g., [31, 61, 103]), Eisenstat’s example
can be computed directly using n/2 + 2 groups by separately evaluating submatrices
A1 and A2 . Clearly, the columns of A1 can be partitioned into two structurally
orthogonal groups with columns of D1 in the ﬁrst and columns of D2 in the second.
For each i, i ≤ n/2, columns i and n/2 + i of A2 are structurally orthogonal, andhence by grouping two such columns together, a structurally orthogonal partition of
A2 consisting of n/2 groups is possible. Hence, the entries of matrix A can be obtained
from n/2 + 2 matrix-vector products. An evaluation of A in this manner corresponds
to using = 2 in the segmented column approach of Hossain and Steihaug.
For the case where n = 6 the discussion above shows that Eisenstat’s example
can be computed directly using ﬁve, instead of six, matrix-vector products. Hossain
and Steihaug [65] have recently shown that it can be computed using four groups
and that four is the optimal number of groups required in a direct method. They
considered the column grouping {a3 , a4 , a5 }, {a2 , a4 , a6 }, {a1 , a5 , a6 }, and {a1 , a2 , a3 }which corresponds to the seed matrix

0
0

1
D=
1

1
0

0
1
0
1
0
1

1
0
0
0
1
1


1
1

1
.
0

0
0

The matrix-vector product AD from which the entries of A can be directly recovered

690

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

is thus



a14

a25


a33

a43
AD = 


a55


a64
a73 + a74 + a75

a14
a22
a36
a42
a56
a62 + a64 + a66
a74

a11
a25
a36
a41
a51 + a55 + a56
a66
a75


a11

a22


a33

a41 + a42 + a43 
.

a51


a62
a73

We note that the EI edge coloring (see the left part in Figure 10.2) in which theedges are partitioned into the four groups {a11 , a41 , a66 , a75 }, {a14 , a25 , a33 , a43 , a55 , a64 },{a22 , a36 , a42 , a56 , a74 }, and {a51 , a62 , a73 } corresponds to the seed matrix of Hossain
and Steihaug.
11. Miscellaneous Topics. Some of the coloring problems addressed in this paper
have been studied outside the context of derivative matrix computation from a purely
graph-theoretic or algorithmic point of view. In this section, we brieﬂy review partof the relevant literature with emphasis on work that is close to the subject of thispaper.
We begin in section 11.1 by discussing examples of eﬀective vertex orderings fora greedy distance-1 coloring heuristic. In section 11.2 we include a relatively detailed
discussion of the coloring number of a graph, a concept also useful in contexts otherthan graph coloring. In section 11.3 we mention a few results on upper and lower
bounds for the chromatic number of a graph. Finally, in section 11.4 we provide an
overview of complexity and approximability results for some of the coloring problems
addressed in this paper. For further information on some of the issues discussed in
sections 11.2–11.4, we refer the reader to Jensen and Toft [71] and Toft [115].
Our discussion in sections 11.1 through 11.3 pertains to distance-1 coloring. In
those sections, for brevity, we write “coloring” instead of “distance-1 coloring.” Insection 11.4, however, we will explicitly state the kind of coloring being discussed.
Unlike previous sections of this paper, the current section assumes familiarity withmore advanced graph-theoretic concepts; we will use several such concepts without
deﬁning them.
11.1. Ordering for Coloring. In 1972, Matula, Marble, and Isaacson [94] analyzed various practically eﬀective greedy coloring heuristics. A decade later Coleman
and Moré [30] used, among others, some of these heuristics to solve a Jacobian estimation problem. Manvel [92], Kubale [85], and very recently Kosowski and Manuszewski[80] have surveyed greedy coloring heuristics.
A key issue in a greedy coloring heuristic is the order in which the vertices are
visited as the latter determines the number of colors used by the heuristic. Here we
discuss a few examples of eﬀective vertex orderings.
For much of the discussion in this section we consider the greedy coloring heuristic
SEQ, short for “sequential,” outlined in Algorithm 11.1.
Notation. We begin by collecting some notation to be used here and in section 11.2.
For a graph G = (V, E) and a given vertex ordering v1 , v2 , . . . , vn , let G[Vi ] denotethe graph induced by the vertex set Vi = {v1 , . . . , vi }. Hence G[Vn ] ≡ G. Let qi
be the number of colors required by SEQ to color G[Vi ]; thus qn is the number ofcolors required by SEQ to color the entire graph G. Let d(v, G[Vi ]) be the number ofvertices adjacent to vertex v in G[Vi ]; thus d(v, G[Vn ]) ≡ d(v) is the degree of vertex
v in G. Recall that ∆(G) and δ(G) denote the maximum and minimum degree inGRAPH COLORING FOR COMPUTING DERIVATIVES

691

Algorithm 11.1. Greedy (sequential) coloring.
procedure SEQ(G = (V, E))
Let v1 , v2 , . . . , vn be a vertex ordering
Assign v1 color 1
for i ← 2 to n do
Assign vi the smallest color not used by any of its neighbors
end forend procedure

G, respectively. When the graph under discussion is clear from the context, we may
simply use ∆ and δ.
Theorems 11.1 and 11.3 below give upper bounds on the number of colors used by
SEQ that reﬂect the importance of the chosen ordering; we will shortly see orderings
that minimize the upper bounds given in the two theorems.
Theorem 11.1. The number of colors used by SEQ satisﬁes the inequalityqn ≤ max {d(vi , G[Vi ])} + 1.
1≤i≤nProof. In the ith step, i > 1, if the algorithm does not introduce a new color tocolor vertex vi , then by deﬁnition of qi , we have qi = qi−1 . Otherwise, qi = qi−1 + 1.
In the latter case, d(vi , G[Vi ]) satisﬁes the inequality qi−1 ≤ d(vi , G[Vi ]) since a newcolor was needed to color vi . The theorem follows by induction on i.
The following result is an immediate consequence of Theorem 11.1.
Corollary 11.2. For any ordering, SEQ uses at most ∆ + 1 colors.
The theorem below gives a weakening of the bound given in Theorem 11.1.
Theorem 11.3. The number of colors used by SEQ satisﬁes the inequalityqn ≤ max min{d(vi ), i − 1} + 1.
1≤i≤nProof. Consider the bound in Theorem 11.1. Clearly d(vi , G[Vi ]) ≤ d(vi , G) ≡d(vi ). Also, since there are exactly i − 1 vertices in G[Vi ] other than vi itself, thestatement d(vi , G[Vi ]) ≤ i − 1 holds—hence the bound in this theorem.
Notice that a nonincreasing ordering of the vertices with respect to their degrees
in the input graph G minimizes the upper bound in Theorem 11.3. This ordering,known as largest ﬁrst (LF) ordering, was suggested by Welsh and Powell [120]. In
general, however, the number of colors used by SEQ using an LF ordering can bemuch larger than the optimal number of colors. An example of a bipartite (hence
2-colorable) graph with n vertices and maximum degree n/2 − 1 for which an LF
ordering could require n/2 colors is given in [30]. The example, as a by-product,shows that the bound of Corollary 11.2 is tight.
An ordering that has been found to usually require fewer colors compared withan LF ordering is a smallest last (SL) ordering. In SL ordering the last vertex in theordering, vn , is chosen to be a vertex with the minimum degree in G. The order ofthe remaining vertices is deﬁned backwards. Assume that the vertices vi+1 , . . . , vn inthe ordering have been determined. Then vertex vi is chosen to be a vertex with theminimum degree in the subgraph induced by the set V − {vi+1 , . . . , vn }. At each step,
while choosing a vertex of minimum degree, ties are broken arbitrarily. Note that
although an SL ordering is computed backwards, vn through v1 , the vertices are thencolored in the order v1 through vn . In this, SL ordering diﬀers from LF ordering.

692

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

Algorithm 11.2. Greedy coloring formulated as a while-loop.
procedure SEQ-W(G = (V, E))U ←V> U is the set of currently uncolored verticesi←1while U = ∅ do
Choose a vertex vi from U using some criterion
Assign vi the smallest color not used by any of its neighbors
U ← U − {vi }i←i+1
end while
end procedure

The SL ordering was proposed by Matula [93] and, as discussed in section 11.2,
its relevance in other contexts was discovered independently by other researchers.
In the context of algorithm SEQ, an SL ordering minimizes the bound given in
Theorem 11.1. In particular, in an SL ordering, for each i, 1 ≤ i ≤ n,d(vi , G[Vi ]) = min d(vj , G[Vi ]).
vj ∈ViAn SL ordering of the vertices in a graph can be obtained in time proportional
to the number of edges in the graph.
In the literature other vertex orderings, such as incidence degree (ID) and saturation degree (SD) ordering, have also been suggested as heuristics for reducing thenumber of colors used by a greedy algorithm. The former ordering is due to Coleman
and Moré [30] and the latter is due to Brélaz [22]. ID and SD orderings are similar in
spirit to LF ordering; they diﬀer in the specialized “degree” measure they use. Unlike LF, SL, and ID orderings, an SD ordering can be computed only as the coloringproceeds. In other words, an SD ordering cannot be precomputed to then color thevertices.
To appreciate the ideas in LF, ID, and SD orderings, consider the scheme SEQW outlined in Algorithm 11.2. At the ith step of SEQ-W that uses LF ordering,
among the current set of uncolored vertices, a vertex having the maximum degreein the input graph is selected and colored. An ID ordering–based coloring improveson this by focusing on the already colored part of the input graph. In particular,
the ith step of SEQ-W that uses ID ordering selects a vertex (from the current setof uncolored vertices) having the maximum incidence degree, the number of alreadycolored neighbors. Note that the incidence degree of vertex vi is precisely d(vi , G[Vi ]).
In a further improvement, the ith step of SEQ-W that uses SD ordering chooses
a vertex having the maximum saturation degree, the number of distinctly coloredneighbors. In each of the three orderings LF, ID, and SD, ties are broken arbitrarily.
The notions saturation degree (sd), incidence degree (id), and usual degree (d) of avertex v satisfy the inequalities sd(v) ≤ id(v) ≤ d(v). Figure 11.1 illustrates this; inthe ﬁgure vertex c has the values two, three, and four for sd, id, and d, respectively.
Notice that a vertex with a larger “degree” is more constrained in the choice of colors
than a vertex with a smaller “degree.” Note also that among the various specialized
degree measures, saturation degree is the measure that most accurately models the
choice of colors available to a vertex. Thus, intuitively one could expect algorithm
SEQ-W to require a nondecreasing number of colors while using SD, ID, and LF
orderings. Experimental results of Coleman and Moré [30] on graphs that arise inGRAPH COLORING FOR COMPUTING DERIVATIVES

b

693

c

d

a

f

e

Fig. 11.1 A partially colored graph illustrating the notions of saturation degree (sd), incidence degree(id), and “usual” degree (d). In this illustration, vertex a is red, vertices b and e are blue,and vertex f is green. Vertex c has sd = 2, id = 3, and d = 4.
Table 11.1 Summary of properties of SL, LF, ID, and SD ordering.

Ordering

ith vertex satisﬁes

SLd(vi , G[Vi ]) = min d(vj , G[Vi ])LFd(vi , G) =ID

d(vi , G[Vi ]) =SD

sd(vi , G[Vi ]) =vj ∈Vimax

vj ∈(V −Vi−1 )d(vj , G)maxd(vj , G[Vi−1 ∪ {vj }])max

sd(vj , G[Vi−1 ∪ {vj }])vj ∈(V −Vi−1 )vj ∈(V −Vi−1 )practice by and large support this intuition. In terms of time complexity, algorithm
SEQ-W using LF or ID ordering can be implemented to run in O(m) time, whereas
an SD ordering–based coloring runs in O(n2 ) time [22, 30, 120]. Here n and m are
the numbers of vertices and edges, respectively.
In Table 11.1 we give a succinct summary of the properties of SL, LF, ID, andSD orderings. The right column of the table shows a statement satisﬁed by the ith
vertex in the respective ordering, for each i, 1 ≤ i ≤ n; in the case where i = 1, theset Vi−1 on the right column corresponds to the empty set. The quantity sd(v, G[Vi ])in the last row of the table shows the saturation degree of the vertex v in the graph
G[Vi ]. Notice that any vertex in G can be chosen to be v1 in the case of ID and SD
orderings.
Coleman and Moré [30] showed that algorithm SEQ-W that uses an ID ordering
requires two colors on bipartite graphs and hence is optimal for this class of graphs.
They also showed that this is not the case for SL ordering by giving an example of abipartite graph on n vertices in which SEQ using an SL ordering could require n/3+1
colors.
It should be pointed out that there exists an ordering of the vertices in which SEQ
would use the optimal number of colors. To see this, consider an optimal coloring
φ : V → {1, 2, . . . , p} of a graph G = (V, E). Let the vertices be ordered such that thesequence {φ(vi )} is nondecreasing, i.e., φ(vi ) ≤ φ(vi+1 ) for 1 ≤ i ≤ n−1. Note that insuch an ordering, vertices of each color class are listed consecutively. It is then easy to
see that if the vertices were to be colored afresh in the order just described, algorithm
SEQ would produce a coloring with an optimal number of colors. The problem while694

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

2
b

3
c

3
b
6
d

1
a

2
c
6
d

4
a

f

e

f

e

4

5

1

5

a

b

Fig. 11.2 Two vertex orderings of a graph with diﬀerent values for maximum back-degree; themaximum back-degree in (a) is three and in (b) is two (which is optimal).

using SEQ, of course, is ﬁnding an ordering among the n! possibilities that would
produce an optimal coloring. The fact that there exists an ordering in which SEQ
would use the optimal number of colors is unlike a number of other problems in
optimization in which a greedy algorithm provided with the “right ordering” might
still fail to yield an optimal solution.
11.2. The Coloring Number. In a given ordering of the vertices of a graph, let
the back-degree of a vertex vi be the number of adjacent vertices to vi that precede vi inthe ordering. Using the notation introduced in section 11.1, the back-degree of vertex
vi is the same as d(vi , G[Vi ]), a quantity that is also referred to as the incidence degreeof vertex vi . For a given graph G = (V, E), let Bπ (G) denote the maximum backdegree in vertex ordering π : V ↔ {1, 2, . . . , n}. Clearly, diﬀerent vertex orderings
may result in diﬀerent maximum back-degrees. As an illustration, Figure 11.2 shows
two diﬀerent vertex orderings in a graph. The maximum back-degree in the ordering
shown in part (a) is three and the corresponding value in part (b) is two.
The quantity minπ {Bπ (G)}, where the minimum is taken over the n! possible
vertex orderings, is of special interest. Let this quantity be denoted by B ∗ (G). For thegraph shown in Figure 11.2, the ordering shown in part (b) is optimal, i.e., B ∗ (G) = 2.
Consider the coloring algorithm SEQ in which the vertex ordering used is suchthat B ∗ (G) is attained. With such an ordering, in each step i of SEQ, vertex vi has at
most B ∗ (G) already colored neighbors and hence can be colored using one value fromthe set {1, 2, . . . , B ∗ (G) + 1}. Thus B ∗ (G) + 1 colors suﬃce to color G. Motivated bythis fact, the quantity B ∗ (G) + 1 is known in the literature as the coloring number of
G and is denoted by col(G). The name coloring number was ﬁrst used by Erdös andHajnal [40]. As we will shortly see, the graph parameter B ∗ (G) is also of interest incontexts other than coloring.
The coloring number may appear to be an intractable graph parameter as theminimum is evaluated over n! possible orderings. However, this was proven to be
not the case independently by Finck and Sachs [42] and Matula [93]. Theorem 11.4
reﬂects this fact. It states that the number of colors required by algorithm SEQ on agraph G with an SL ordering is the coloring number of G. Recall that an SL ordering
can be computed in linear time in the number of edges in a graph.
Theorem 11.4. Every graph G satisﬁescol(G) = BSL (G) + 1.

GRAPH COLORING FOR COMPUTING DERIVATIVES

695

Proof. By deﬁnition, col(G) = B ∗ (G) + 1. Clearly B ∗ (G) ≤ BSL (G). Hencecol(G) ≤ BSL (G)+1. We now show that the inequality holds in the opposite direction
as well, thus implying equality.
Consider an SL ordering S = v1 , v2 , . . . , vn of the vertices of G. Clearly, for every
i < n, col(G) ≡ col(G[Vn ]) ≥ col(G[Vi ]); furthermore, col(G[Vi ]) ≥ δ(G[Vi ]) + 1,
since the back-degree of vertex vi in the subsequence v1 , v2 , . . . , vi of S is its ordinary
degree in G[Vi ], which in turn is equal to δ(G[Vi ]) since S is an SL ordering. Hencecol(G) ≥ maxi {δ(G[Vi ])} + 1. The quantity on the right side of the inequality in thepreceding statement is precisely BSL (G) + 1. Thus we have col(G) ≥ BSL (G) + 1,
which completes the proof.
Matula [93] and Szekeres and Wilf [113], independently, have related the coloring number of a graph to yet another seemingly intractable graph parameter. The
parameter we are referring to is the maximum minimum degree in an induced subgraph of a graph, where the maximum is taken over all possible induced subgraphs.
In a given graph G = (V, E), since there are 2n possible subsets of V , computing the
maximum minimum degree in an induced subgraph of G appears to be intractable.
Theorem 11.5 shows this is not the case. Let H ⊆ G denote a nonempty induced
subgraph of G.
Theorem 11.5. Every graph G satisﬁescol(G) = max {δ(H)} + 1.
H⊆G

Proof. Let k = maxH⊆G {δ(H)}. We claim that the graph G = (V, E) itself has avertex of degree at most k, for otherwise the minimum degree in G, δ(G), is strictly
greater than k, contradicting the deﬁnition of k. Let vn be a vertex with degree ≤ kin the graph G and let Hn−1 = G[V − {vn }]. By assumption, Hn−1 has a vertex
of degree at most k. Let vn−1 be one of them and let Hn−2 = G[V − {vn , vn−1 }].
Continue in this manner to obtain the ordering v1 , v2 , . . . , vn of the vertices of G.
It is then clear that in such an ordering the maximum back-degree of a vertex is at
most k. Hence col(G) ≤ k + 1.
We now prove that the inequality holds in the other direction as well. For any
H ⊆ G, clearly col(G) ≥ col(H) and col(H) ≥ δ(H) + 1, since the back-degree ofthe last vertex in any ordering of the vertices of H is just its ordinary degree in H,
which is at least δ(H). Thus col(G) ≥ k + 1. The theorem follows since col(G) isshown to be both less than or equal to and greater than or equal to the quantity
maxH⊆G {δ(H)} + 1.
Clearly, the coloring number of a graph is a tighter upper bound on its chromatic
number than the bound ∆ + 1; i.e., for any graph G the inequalities χ(G) ≤ col(G) ≤∆(G) + 1 hold. Moreover, col(G) = ∆(G) + 1 if and only if G has a ∆(G)-regularconnected component.
A graph G for which col(G) ≤ k + 1 has also been called k-degenerate by Lick andWhite [88]. (A graph is k-degenerate if every induced subgraph has minimum degree
at most k.) A related concept has been studied in the social and biological network
literature as the k-core of a graph [21, 105]. The k-core of a graph is a maximal
induced subgraph in which every vertex has at least k neighbors in the subgraph. A
k-core of a graph can be computed by a linear-time algorithm that visits the vertices
in the reverse of an SL ordering. A paper by Ramadan, Tarafdar, and Pothen [109]extends the deﬁnition of the k-core to a hypergraph and discusses the k-cores ofprotein interaction graphs and protein complex hypergraphs. The coloring number ofa graph is also closely related to the arboricity of the graph [37].

696

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

11.3. Bounds for the Chromatic Number. Since computing the chromatic number χ(G) of a general graph G is NP-hard, we consider bounding it with more tractablegraph parameters.
Upper Bounds. As we have already seen, the coloring number col(G) is a lineartime computable upper bound for the chromatic number χ(G) of a graph G. We havealso seen that ∆(G) + 1 is another easy to compute (but weaker) upper bound for
χ(G). If G is a complete graph or an odd cycle, then in fact χ(G) = ∆(G) + 1. In
all other cases, Brooks’s theorem tells us that the general bound ∆(G) + 1 can be
improved by 1.
Theorem 11.6 (Brooks’s theorem). Let G be a connected graph. If G is neither
complete nor an odd cycle, then χ(G) ≤ ∆(G).
For a proof of Brooks’s theorem, see, e.g., Diestel [37].
Lower Bounds. The size of the largest clique in G, ω(G), is one obvious lower
bound for χ(G). Another, less obvious, lower bound comes from the fact that agraph coloring is a vertex partition into independent sets. If G does not contain an
independent set of size s+1, then in every coloring of G at most s vertices get the samecolor. Hence, χ(G) ≥ |V |/α(G), where α(G) is the size of the largest independent set
in G = (V, E). Combining the two, one gets the following improved lower bound:χ(G) ≥ max{ω(G), |V |/α(G)}.
For many graphs the bound just given is a weak one. Moreover, both ω(G) andα(G) are NP-hard to determine. By considering a more general class of subgraphs
than cliques, Hajós [57] obtained a further improved lower bound (whose computation
is not known to be polynomial). In particular, Hajós showed that χ(G) ≥ k if andonly if G has a k-constructible subgraph. For further discussion on Hajós’s and otherlower bounds, see, e.g., Diestel [37] or Toft [115].
11.4. Theoretical Results on Coloring Problems. In this last subsection, we
provide an overview of various results in the graph theory literature on some of thecoloring problems addressed in this paper. Our main objective here is to point theinterested reader to relevant references.
Each of the NP-hardness results stated in the next paragraph has been mentioned
elsewhere in this paper; here the results are presented together for convenience.
Lin and Skiena [89] proved that the distance-k graph coloring problem is NP-hardfor every ﬁxed integer k ≥ 1. The proof relies on the equivalence between distance-kcoloring of G and distance-1 coloring of Gk . Coleman and Moré [31] showed that thestar coloring problem is NP-hard even if the graph is bipartite. Coleman and Verma[32] showed that the problem of ﬁnding a star bicoloring using the fewest colors is
also NP-hard. Further, the acyclic coloring problem was proven to be NP-hard by
Coleman and Cai [27].
For a general graph on n vertices, the ﬁrst approximation algorithm for distance-1
coloring, with an approximation ratio of O(n log1 n ), was obtained by Johnson [72]. This2
3
was later improved to O(n( logloglogn n ) ) by Wigderson [121] and to O(n( logloglogn n ) ) by2

log n)),
Berger and Rompel [14]. The current best known approximation ratio is O(n (log(log n)3a result due to Halldórsson [59]. Based on Wigderson’s algorithm, Karger, Motwani,and Sudan [75] considered approximate graph coloring using semideﬁnite programming techniques. On the negative side, Bellare, Goldreich, and Sudan [12] showedthat distance-1 coloring is not approximable within O(n1/7− ) for any > 0, unless P= NP.

GRAPH COLORING FOR COMPUTING DERIVATIVES

697

A planar graph can be distance-1 colored using four colors. This fact, best known
as the four color theorem, is perhaps one of the most famous results in graph theory.
After challenging mathematicians for nearly a century, the theorem was proved byAppel and Haken in 1976 [7, 10]. However, due to its length, its extensive use ofcomputer veriﬁcation, and its omission of details, the proof of Appel and Haken has
been subject to criticism. In response, the authors published a revised version of theproof a decade later [8, 9]. Jensen and Toft [71] cited a 1994 manuscript of Robertson,
Sanders, Seymour, and Thomas that presents a new, highly simpliﬁed proof of thefour color theorem. For further information on the four color theorem, see the recent
books by Fritsch and Fritsch [43] and Wilson [122].
Since every planar graph has a vertex of degree at most ﬁve, one can give a
linear-time algorithm for ﬁnding a distance-1 coloring of a planar graph that uses atmost six colors: run algorithm SEQ with an SL ordering. A linear-time algorithm
that ﬁnds a distance-1 coloring of a planar graph using ﬁve colors was presented by
Matula, Shiloach, and Tarjan [95, 96] and independently by Chiba, Nishizeki, andSaito [25]. Also, a linear-time algorithm that uses ﬁve colors is implicit in a recent
proof of the theorem “every planar graph is ﬁve colorable” given by Thomassen [114].
Deciding whether a planar graph of maximum degree four can be distance-1 coloredusing three colors is NP-complete [71].
Wegner [119] studied distance-k (k ≥ 2) chromatic numbers of planar graphs
having maximum degree at most three and proved that for k = 2 eight colors suﬃce,
whereas seven is the best possible. Recently, Agnarsson and Halldórsson [3] studied
distance-k coloring of planar graphs. They showed that a planar graph having a sufﬁciently large maximum degree ∆ can be distance-2 colored using at most 9∆/5
colors. The authors also showed that for a ﬁxed integer k ≥ 1, the distance-k chromatic number of a planar graph is O(∆k/2 ). Agnarsson, Greenlaw, and Halldórsson[2] studied distance-k coloring of chordal graphs and showed that for even integers k,
the problem is not approximable within a factor of |V |1/2− for any > 0.
Acyclic colorings were ﬁrst deﬁned and studied by Grünbaum [56]. Grünbaum
asked whether the acyclic chromatic number χa (G) of a graph G satisﬁes the inequality χa (G) ≤ ∆ + 1. A negative answer to this question was given by Erdös, whoproved probabilistically the existence of graphs where χa (G) ≥ ∆4/3− [71]. Alon,
McDiarmid, and Reed [6] showed (again by probabilistic methods) that there exists a
constant c1 > 0 such that χa (G) ≤ c1 ·∆4/3 . The same authors also showed that there−1/3
is a constant c2 > 0 such that there exist graphs with χa (G) > c2 · ∆4/3 · (log ∆)for inﬁnitely many ∆. There are a few studies concerning χa (G) for speciﬁc values of∆. It is clear that χa (G) ≤ 2 if and only if G is a forest. Grünbaum [56] proved that
χa (G) ≤ 4 for ∆ = 3, and χa (G) ≤ 6 for ∆ = 4. The latter bound was reduced fromsix to ﬁve by Burstein [23]. Kostochka [82] proved that it is NP-complete to decide
whether a given graph G satisﬁes χa (G) ≤ 3.
Grünbaum [56] also proved that every planar graph admits an acyclic coloring
using at most nine colors. After a series of improvements [4, 81, 98] the number was
ﬁnally reduced to ﬁve by Borodin [19]. Borodin’s proof is similar to the proof of
the four color theorem in that there are 450 reducible conﬁgurations that need to bechecked.
Star coloring has been recently studied by Albertson et al. [5] who show that every
acyclic coloring that uses q colors can be reﬁned to a star coloring with at most 2q 2 −q
colors. They also prove that planar graphs have star colorings with at most 20 colors
and exhibit an example in which 10 colors are required. This is an improvement on
the upper bounds 80 and 30, obtained earlier by Fertin, Raspaud, and Reed [41] and698

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

Neśetŕil and Ossan de Mendez [102], respectively. The problem of deciding whether
a planar bipartite graph has a star coloring with at most three colors was shown to
be NP-complete in [5].
12. Conclusion.
12.1. Summary. We have studied the eﬃcient computation of sparse Jacobianand Hessian matrices using ﬁnite diﬀerence and automatic diﬀerentiation techniques.
In this context, we considered a variety of matrix partitioning problems that can beclassiﬁed along four axes. The ﬁrst axis corresponds to whether the matrix to be computed is nonsymmetric (Jacobian) or symmetric (Hessian). The second axis relates towhether the matrix partition used is unidirectional (involving either columns or rows)
or bidirectional (involving both columns and rows). The third axis corresponds towhether the matrix entries are determined directly (by solving a diagonal system of
equations) or via substitution (by solving a triangular system of equations). The lastaxis shows whether all the entries of the matrix need to be computed or just a subsetof them suﬃces for preconditioning purposes. Out of the sixteen possibilities this scenario suggests, ten cases correspond to meaningful problems within our framework,and eight of them are addressed in this paper.
Building upon several pioneering papers that have been published in the lasttwenty years, we have developed a unifying graph-theoretic framework for studying
these matrix partitioning problems.
Our problem formulations and algorithms are based on robust and space-eﬃcientgraph representations of the sparsity structure of matrices: The Jacobian is represented by its bipartite graph, and the Hessian by its adjacency graph.
We showed that the distance-2 graph coloring problem is a generic model forthe various matrix partitioning problems. Our approach provides fresh insight into
the matrix partitioning problems, leading to several simple and eﬀective heuristic
algorithms. Among others, we developed fast star coloring heuristics, sketched a
new heuristic algorithm for the acyclic coloring problem, and suggested a scheme forsolving bicoloring problems.
We have also shown a hypergraph coloring formulation for the unidirectional
Jacobian computation problem. This alternative formulation may be helpful in partial
matrix computation for preconditioning purposes.
In the case of unidirectional Jacobian computation via a direct method, we showed
experimentally the advantages oﬀered by the new partial distance-2 coloring formulation as compared to the previously known distance-1 coloring formulation. In computing the Hessian using a direct method, we demonstrated a time/quality trade-oﬀ
between two star coloring algorithms.
In general, most of the algorithms in the literature for solving the coloring problems considered in this paper rely on ﬁrst transforming the input graph G = (V, E)
to some denser graph G = (V, E  ), E  ⊇ E, such that a distance-1 coloring of G
is equivalent to the particular coloring problem on G [27, 31, 32]. In contrast, thealgorithms proposed in this paper solve the particular coloring problem directly on
G. The main advantages oﬀered by our approach include: smaller storage space requirement, the possibility to mix-and-match methods, and ease of developing more
eﬃcient algorithms and ﬂexible software.
Our primary emphasis has been on computational methods that treat a column
or a row of a matrix as an indivisible unit and that partition the columns or rows (orboth). However, we have brieﬂy reviewed methods that do not ﬁt this framework.
We formulated one of these methods, element isolation, as a specialized edge coloringin the bipartite graph of the Jacobian.

GRAPH COLORING FOR COMPUTING DERIVATIVES

699

We have also included a detailed discussion of eﬀective vertex orderings for greedy
coloring, a discussion of the coloring number of a graph and its signiﬁcance in various
contexts, and a brief review of relevant theoretical results from the graph coloring
literature.
12.2. Recent Development. During the time this paper was under review or inpress, we have made progress on two fronts. First, we have implemented the acyclic
coloring algorithm sketched in section 6.1.3 and developed a new algorithm for thestar coloring problem; both algorithms are based on the paradigm of maintaining twocolored subgraphs. A paper that discusses these algorithms, and their implementationand performance has been submitted elsewhere [48]. Our implementations of the
new acyclic and star coloring algorithms are currently being incorporated into theautomatic diﬀerentiation software ADOL-C [118]. Second, we have implemented thevarious vertex orderings discussed in section 11.1, each tailored for the distance-k(k = 1, 2) coloring problem. In a future paper, we plan to report experimental results
demonstrating the impact of the various ordering techniques when sequential coloringalgorithms are applied to large-scale problems. All software we develop will be madepublicly available.
12.3. Further Work. This work can be followed up along several directions,some of which we list below. We hope that these problems will attract the attentionof our readers.
1. The new algorithms for star bicoloring (section 5.4) and acyclic bicoloring(section 6.2.3) suggested in this paper need to be implemented.
2. New algorithms need to be designed and implemented for the newly considered partial matrix computation problems (section 8).
3. One of the motivations for the current study has been the need to develop
parallel algorithms for solving partitioning problems in large-scale optimization contexts. In previous work [45, 47] we developed shared-memory parallel
algorithms for the distance-2 and star coloring problems. Some of the ideas
in these algorithms have recently been extended to design and implement
distributed memory parallel algorithms for the distance-1 and the distance-2
coloring problems [18, 20]. Based on these early experiences and the current study, scalable distributed-memory parallel algorithms for the various
coloring problems need to be designed and implemented.
4. In the case of bidirectional partitioning problems, based on the known relationship to graph bicoloring, we have argued that ﬁnding a “small”-size vertex
cover in a preprocessing step contributes to making the overall computation
more eﬃcient (section 5.4). Loosely speaking, the desired vertex cover has
to favor high-degree vertices from both vertex sets of the bipartite graph. A
precise characterization of the “optimum” vertex cover required is an issue
that we are currently studying.
5. The maximum number of nonzeros in a row of a Jacobian is a lower boundon the number of colors required to distance-2 color the column vertices of
the associated bipartite graph Gb = (V1 , V2 , E). A tighter lower bound is thesize of a maximal clique in the graph G2b [V2 ] (section 3.5.2). How can one
ﬁnd such a clique and its size quickly?
6. What nontrivial lower bounds can be given for star coloring (section 4)? forstar bicoloring (section 5)? for acyclic coloring (section 6.1)? for acyclic
bicoloring (section 6.2)?
7. Are there special graph classes for which these problems are tractable or havegood upper bounds?

700

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN

8. This paper dealt with coloring formulations of Jacobian and Hessian computation problems where the sparsity pattern is not necessarily regular. In
cases where the pattern is regular, such as Jacobians or Hessians that arise
in discretization of structured grids, the resulting coloring problems can besolved optimally in polynomial time (see, e.g., the work of Goldfarb andToint [52]). A closer look at coloring problems for such structured graphs isa worthwhile issue.
9. For the variant of edge coloring discussed in section 10.3, the maximum degreein the set of row vertices of the bipartite graph is a lower bound for the fewest
colors needed. What upper bound can be given? What exact/approximation
algorithm can be designed?
10. Star coloring a planar graph requires at most twenty colors, and an example
is known for which ten colors are needed (section 11.4). Reducing the gap
between the lower bound and the upper bound is an open problem.
Graph coloring for eﬃciently computing Jacobians and Hessians represents oneof the earliest applications of combinatorial methods in numerical optimization. Thisis but one of the problems in the emerging research area of combinatorial scientiﬁccomputing, an area where combinatorial models and algorithms are used to solve
problems in scientiﬁc computing. Our broader hope is that this article provides animpetus for the formulation and solution of many more combinatorial problems inscientiﬁc computing.
Acknowledgments. This work was begun when Assefaw Gebremedhin visited
Alex Pothen at Old Dominion University for a semester during his Ph.D. study at the
University of Bergen, Norway; Fredrik Manne was Assefaw’s advisor. Paul Hovland
and Rob Bisseling, who served on Assefaw’s thesis examination committee, read earlier versions of some of the sections of this paper and suggested several improvements
and new problems to consider. Rob Bisseling has also read a recent version of thispaper and shared his constructive remarks with us. Trond Steihaug at the Universityof Bergen has given us many helpful comments. Tom Coleman encouraged our workon graph coloring for optimization, even when others had felt that “it has all been
done already.” Paul Hovland and Lois McInnes have been enthusiastic about using
this work in their software for AD and optimization. Our work on partial estimation
problems were stimulated by conversations with Bill Spotz at Sandia National Labs.
Bruce Hendrickson, Erik Boman, Robert Preis (currently at the University of Paderborn, Germany), and other members of the Discrete Algorithms Group at Sandia
National Labs have discussed these results with us during Alex Pothen’s sabbatical
visits to Sandia. Dinesh Kaushik gave us some insight on the runtime of a ﬁnite
diﬀerence code for computing the numerical values of the entries of a large, sparse
Jacobian matrix. Randy LeVeque and Michael Overton have read an earlier version ofthis paper and given us helpful feedback and encouragement. The anonymous referees
have provided helpful comments on presentation. This paper has been shaped andimproved by our interactions with these supportive colleagues, and our thanks go to
all of them.
REFERENCES[1] A. Aggarwal, A. Bar-Noy, D. Coppersmith, R. Ramaswami, B. Schieber, and M. Sudan,Eﬃcient routing and scheduling algorithms for optical networks, in Proceedings of the
5th Annual ACM-SIAM Symposium on Discrete Algorithms, SIAM, Philadelphia, 1994,
pp. 412–423.

GRAPH COLORING FOR COMPUTING DERIVATIVES

701[2] G. Agnarsson, R. Greenlaw, and M. M. Halldórsson, On powers of chordal graphs andtheir colorings, Congr. Numer., 100 (2000), pp. 41–65.
[3] G. Agnarsson and M. M. Halldórsson, Coloring powers of planar graphs, SIAM J. Discrete
Math., 16 (2003), pp. 651–662.
[4] M. O. Albertson and D. M. Berman, The acyclic chromatic number, Congr. Numer., 17(1976), pp. 51–69.
[5] M. O. Albertson, G. G. Chappell, H. A. Kierstead, A. Kündgen, and R. Ramamurthi,Coloring with no 2-colored P4’s, Electron. J. Combin., 11 (2004), article R26.
[6] N. Alon, C. McDiarmid, and B. Reed, Acyclic coloring of graphs, Random Structures
Algorithms, 2 (1991), pp. 277–288.
[7] K. Appel and W. Haken, Every planar map is four colorable, Part I: Discharging, Illinois
J. Math., 21 (1977), pp. 429–490.
[8] K. Appel and W. Haken, The four color proof suﬃces, Math. Intell., 8 (1986), pp. 10–20.
[9] K. Appel and W. Haken, Every Planar Map Is Four Colorable, AMS, Providence, RI, 1989.
[10] K. Appel, W. Haken, and J. Koch, Every planar map is four colorable, Part II: Reducibility,
Illinois J. Math., 21 (1977), pp. 491–567.
[11] O. Axelsson and U. Nävert, On a graphical package for nonlinear partial diﬀerential equation problems, in Proceedings of the IFIP Congress 77, B. Gilchrist, ed., Information
Processing, North-Holland, Amsterdam, 1977, pp. 103–108.
[12] M. Bellare, O. Goldreich, and M. Sudan, Free bits, PCPs and nonapproximability—
towards tight results, SIAM J. Comput., 27 (1998), pp. 804–915.
[13] C. Berge, Hypergraphs, North-Holland, Amsterdam, 1989.
[14] B. Berger and J. Rompel, A better performance guarantee for approximate graph coloring,
Algorithmica, 5 (1990), pp. 459–466.
[15] J. Bialogrodzki, Path coloring and routing in graphs, in Graph Colorings, M. Kubale, ed.,Contemp. Math. 352, AMS, Providence, RI, 2004, pp. 139–152.
[16] C. Bischof, A. Carle, P. Khademi, and A. Mauer, ADIFOR 2.0: Automatic diﬀerentiationof Fortran 77 programs, IEEE Comput. Sci. Engrg., 3 (1996), pp. 18–32.
[17] C. Bischof, L. Roh, and A. Mauer, ADIC: An extensible automatic diﬀerentiation tool for
ANSI-C, Software—Practice and Experience, 27 (1997), pp. 1427–1456.
[18] E. Boman, D. Bozdağ, U. Catalyurek, A. H. Gebremedhin, and F. Manne, A scalable
parallel graph coloring algorithm for distributed memory computers, in Proceedings ofEuro-Par 2005 Parallel Processing, J. Cunha and P. Medeiros, eds., Lecture Notes in
Comput. Sci. 3648, Springer-Verlag, Berlin, 2005, pp. 241–251.
[19] O. V. Borodin, On acyclic colorings of planar graphs, Discrete Math., 25 (1979), pp. 211–236.
¨ üner,[20] D. Bozdağ, U. Catalyurek, A. H. Gebremedhin, F. Manne, E. Boman, and F. Ozg
A parallel distance-2 graph coloring algorithm for distributed memory computers, in Proceedings of HPCC-05, the 2005 International Conference on High Performance Computing and Communications, L.T. Yang et al., eds., Lecture Notes in Comput. Sci. 3726,
Springer-Verlag, Berlin, 2005, pp. 796–806.
[21] B. J. Breitkreutz, C. Stark, and M. Tyers, Osprey: A network visualization system,
Genome Bio., 4 (2003), article R22.
[22] D. Brélaz, New methods to color the vertices of a graph, Comm. ACM, 22 (1979), pp. 251–256.
[23] M. I. Burstein, Every 4-valent graph has an acyclic 5-coloring, Soobsc. Akad. Nauk Gruzin.
SSR, 93 (1979), pp. 21–24 (in Russian, with Georgian and English summaries).
[24] G. J. Chaitin, M. Auslander, A. K. Chandra, J. Cocke, M. E. Hopkins, and P. Markstein, Register allocation via coloring, Comput. Lang., 6 (1981), pp. 47–57.
[25] N. Chiba, T. Nishizeki, and N. Saito, A linear 5-coloring algorithm for planar graphs, J.
Algorithms, 2 (1981), pp. 317–327.
[26] T. F. Coleman, Large Sparse Numerical Optimization, Lecture Notes in Comput. Sci. 165,
Springer-Verlag, New York, 1984.
[27] T. F. Coleman and J.-Y. Cai, The cyclic coloring problem and estimation of sparse Hessian
matrices, SIAM J. Algebraic Discrete Methods, 7 (1986), pp. 221–235.
[28] T. F. Coleman, B. Garbow, and J. J. Moré, Software for estimating sparse Jacobianmatrices, ACM Trans. Math. Software, 10 (1984), pp. 329–347.
[29] T. F. Coleman, B. Garbow, and J. J. Moré, Software for estimating sparse Hessian matrices, ACM Trans. Math. Software, 11 (1985), pp. 363–377.
[30] T. F. Coleman and J. J. Moré, Estimation of sparse Jacobian matrices and graph coloring
problems, SIAM J. Numer. Anal., 20 (1983), pp. 187–209.
[31] T. F. Coleman and J. J. Moré, Estimation of sparse Hessian matrices and graph coloring
problems, Math. Program., 28 (1984), pp. 243–270.

702

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN[32] T. F. Coleman and A. Verma, The eﬃcient computation of sparse Jacobian matrices using
automatic diﬀerentiation, SIAM J. Sci. Comput., 19 (1998), pp. 1210–1233.
[33] G. Corliss, C. Faure, A. Griewank, L. Hascoët, and U. Naumann, Automatic Diﬀerentiation of Algorithms: From Simulation to Optimization, Springer-Verlag, New York,2002.
[34] J. Culberson, Graph coloring page, http://www.cs.ualberta.ca/∼joe/Coloring/index.html.
[35] A. R. Curtis, M. J. D. Powell, and J. K. Reid, On the estimation of sparse Jacobianmatrices, J. Inst. Math. Appl., 13 (1974), pp. 117–119.
[36] T. Davis, University of Florida Sparse Matrix Collection, NA Digest, 97 (1997); http://www.
cise.uﬂ.edu/research/sparse/matrices.
[37] R. Diestel, Graph Theory, 2nd ed., Springer-Verlag, New York, 2000.
[38] P. Duchet, Hypergraphs, in Handbook of Combinatorics, Volume I, R. L. Graham, M. Grötschel, and L. Lovász, eds., MIT Press, Cambridge, MA, 1996, pp. 381–432.
[39] E. Elmroth, F. Gustavson, I. Jonsson, and B. Kågström, Recursive blocked algorithms
and hybrid data structures for dense matrix library software, SIAM Rev., 46 (2004),
pp. 3–45.
[40] P. Erdös and A. Hajnal, On chromatic number of graphs and set systems, Acta. Math.
Acad. Sci. Hungar., 17 (1966), pp. 61–99.
[41] G. Fertin, A. Raspaud, and B. Reed, On star coloring of graphs, in Graph-Theoretic
Concepts in Computer Science, 27th International Workshop, Lecture Notes in Comput.
Sci. 2204, Springer-Verlag, Berlin, 2001, pp. 140–153.
[42] H. J. Finck and H. Sachs, Über eine von H.S. Wilf angegebene Schranke für die chromatische
Zahl endlicher Graphen, Math. Nachr., 39 (1969), pp. 373–386.
[43] R. Fritsch and G. Fritsch, The Four Color Theorem: History, Topological Foundations,and Idea of Proof, Springer-Verlag, New York, 1998.
[44] A. Gamst, Some lower bounds for a class of frequency assignment problems, IEEE Trans.
Vehicular Tech., 35 (1986), pp. 8–14.
[45] A. H. Gebremedhin, Practical Parallel Algorithms for Graph Coloring Problems in Numerical Optimization, Ph.D. thesis, University of Bergen, Norway, 2003.
[46] A. H. Gebremedhin and F. Manne, Scalable parallel graph coloring algorithms, Concurrency:
Practice and Experience, 12 (2000), pp. 1131–1146.
[47] A. H. Gebremedhin, F. Manne, and A. Pothen, Parallel distance-k coloring algorithms fornumerical optimization, in Euro-Par 2002 Parallel Processing, B. Monien and R. Feldmann, eds., Lecture Notes in Comput. Sci. 2400, Springer-Verlag, Berlin, 2002, pp. 912–
921.
[48] A. H. Gebremedhin, A. Tarafdar, F. Manne, and A. Pothen, New Acyclic and Star Coloring Algorithms with Application to Computing Hessians, Tech. Report, Old Dominion
University, Norfolk, VA, 2005.
[49] U. Geitner, J. Utke, and A. Griewank, Automatic computation of sparse Jacobians byapplying the method of Newsam and Ramsdell, in Computational Diﬀerentiation: Techniques, Applications, and Tools, M. Berz, C. Bischof, G. Corliss, and A. Griewank, eds.,
SIAM, Philadelphia, 1996, pp. 161–172.
[50] R. Giering and T. Taminski, Recipes for Adjoint Code Construction, Tech. Report 212,
Max-Planck Institut für Meteorologie, Hamburg, 1996.
[51] R. K. Gjertsen Jr., M. T. Jones, and P. E. Plassmann, Parallel heuristics for improved,balanced graph colorings, J. Parallel Distrib. Comput., 37 (1996), pp. 171–186.
[52] D. Goldfarb and P. L. Toint, Optimal estimation of Jacobian and Hessian matrices thatarise in ﬁnite diﬀerence calculations, Math. Comp., 43 (1984), pp. 69–88.
[53] A. Griewank, Evaluating Derivatives: Principles and Techniques of Algorithmic Diﬀerentiation, SIAM, Philadelphia, 2000.
[54] A. Griewank and G. F. Corliss, eds., Automatic Diﬀerentiation of Algorithms: Theory,
Implementation, and Application, SIAM, Philadelphia, 1991.
[55] A. Griewank, D. Juedes, and J. Utke, ADOL-C: A package for the automatic diﬀerentiation of algorithms written in C/C++, ACM Trans. Math. Software, 22 (1996), pp. 131–
167.
[56] B. Grünbaum, Acyclic colorings of planar graphs, Israel J. Math., 14 (1973), pp. 390–408.
[57] G. Hajós, Über eine Konstruktion nicht n-färbbarer Graphen, Wiss. Z. Martin-Luther-Univ.
Halle-Wittenberg Math.-Nat. Reihe, 10 (1961), pp. 116–117.
[58] J. Hall, J. Hartline, A. R. Karlin, J. Saia, and J. Wilkes, On algorithms for eﬃcient
data migration, in Proceedings of the 12th Annual ACM-SIAM Symposium on Discrete
Algorithms, SIAM, Philadelphia, 2001, pp. 620–629.
[59] M. M. Halldórsson, A still better performance guarantee for approximate graph coloring,
Inform. Process. Lett., 45 (1993), pp. 19–23.

GRAPH COLORING FOR COMPUTING DERIVATIVES

703[60] J. L. Hennessy and D. A. Patterson, Computer Architecture: A Quantitative Approach,
3rd ed., Morgan-Kaufmann, San Francisco, 2002.
[61] S. Hossain, On the Computation of Sparse Jacobian Matrices and Newton Steps, Ph.D.
thesis, Tech. Report 146, Department of Informatics, University of Bergen, Norway, 1998.
[62] S. Hossain and T. Steihaug, Computing a sparse Jacobian matrix by rows and columns,
Optim. Methods Softw., 10 (1998), pp. 33–48.
[63] S. Hossain and T. Steihaug, Reducing the number of AD passes for computing a sparse
Jacobian matrix, in Automatic Diﬀerentiation of Algorithms: From Simulation to Optimization, G. Corliss, C. Faure, A. Griewank, L. Hascoët, and U. Naumann, eds., SpringerVerlag, New York, 2002, pp. 263–270.
[64] S. Hossain and T. Steihaug, Sparsity issues in computation of Jacobian matrices, in Proceedings of the 2002 International Symposium on Symbolic and Algebraic Computation,
T. Mora, ed., ACM, New York, 2002, pp. 123–130.
[65] S. Hossain and T. Steihaug, Optimal Direct Determination of Sparse Jacobian Matrices,
Tech. Report 254, Department of Informatics, University of Bergen, Norway, 2003.
[66] P. D. Hovland, Personal communication, 2004.
[67] P. D. Hovland, Combinatorial problems in automatic diﬀerentiation, presented at the SIAM
Workshop on Combinatorial Scientiﬁc Computing, 2004.
[68] P. D. Hovland, Automatic Diﬀerentiation of Parallel Programs, Ph.D. thesis, University of
Illinois at Urbana-Champaign, Urbana, IL, 1997.
[69] P. D. Hovland and L. C. McInnes, Parallel simulation of compressible ﬂow using automatic
diﬀerentiation and PETSc, Parallel Comput., 27 (2001), pp. 503–519.
[70] D. Hysom and A. Pothen, A scalable parallel algorithm for incomplete factor preconditioning, SIAM J. Sci. Comput., 22 (2001), pp. 2194–2215.
[71] T. Jensen and B. Toft, Graph Coloring Problems, Wiley-Interscience, New York, 1995.
[72] D. S. Johnson, Worst case behaviour of graph coloring algorithms, Congr. Numer., 10 (1974),
pp. 513–527.
[73] M. T. Jones and P. E. Plassmann, A parallel graph coloring heuristic, SIAM J. Sci. Comput.,
14 (1993), pp. 654–669.
[74] M. T. Jones and P. E. Plassmann, Scalable iterative solution of sparse linear systems,
Parallel Comput., 20 (1994), pp. 753–773.
[75] D. Karger, R. Motwani, and M. Sudan, Approximate graph coloring by semideﬁnite programming, J. ACM, 45 (1998), pp. 246–265.
[76] D. K. Kaushik, Personal communication, 2004.
[77] D. K. Kaushik, D. E. Keyes, and B. F. Smith, On the interaction of architecture andalgorithm in the domain-based parallelization of an unstructured grid incompressible ﬂow
code, in Proceedings of the 10th International Conference on Domain Decomposition
Methods, AMS, Providence, RI, 1998, pp. 311–319.
[78] D. E. Keyes, P. D. Hovland, L. C. McInnes, and W. Samyono, Using automatic differentiation for second-order matrix-free methods in PDE-constrained optimization, in
Automatic Diﬀerentiation of Algorithms: From Simulation to Optimization, G. Corliss,
C. Faure, A. Griewank, L. Hascoët, and U. Naumann, eds., Springer-Verlag, New York,
2001, pp. 35–50.
[79] D. A. Knoll and D. E. Keyes, Jacobian-free Newton-Krylov methods: A survey of approaches and applications, J. Comput. Phys., 193 (2004), pp. 357–397.
[80] A. V. Kosowski and K. Manuszewski, Classical coloring of graphs, in Graph Colorings,
M. Kubale, ed., Contemp. Math. 352, AMS, Providence, RI, 2004, pp. 1–19.
[81] A. V. Kostochka, Acyclic 6-coloring of planar graphs, Metody Diskret. Analiz., 28 (1976),
pp. 40–56 (in Russian).
[82] A. V. Kostochka, Upper Bounds of Chromatic Functions of Graphs, Ph.D. thesis, Novosibirsk, 1978 (in Russian).
[83] M. Krivelevich and B. Sudakov, The chromatic number of random hypergraphs, Random
Structures Algorithms, 12 (1998), pp. 381–483.
[84] S. O. Krumke, M. V. Marathe, and S. S. Ravi, Models and approximation algorithms forchannel assignment in radio networks, Wireless Networks, 7 (2001), pp. 567–574.
[85] M. Kubale, Graph colouring, in Encyclopedia of Microcomputers, Vol. 8, A. Kent and
J. Williams, eds., Marcel Dekker, New York, 1991, pp. 47–69.
[86] M. Kubale, Graph Colorings, AMS, Providence, RI, 2004.
[87] S. L. Lee and P. D. Hovland, Sensitivity analysis using parallel ODE solvers and automatic
diﬀerentiation in C: SensPVODE and ADIC, in Automatic Diﬀerentiation of Algorithms:
From Simulation to Optimization, G. Corliss, C. Faure, A. Griewank, L. Hascoët, andU. Naumann, eds., Springer-Verlag, New York, 2001, Chap. 26, pp. 223–229.

704

A. H. GEBREMEDHIN, F. MANNE, AND A. POTHEN[88] D. R. Lick and A. T. White, k-degenerate graphs, Canad. J. Math., 22 (1970), pp. 1082–
1096.
[89] Y. L. Lin and S. S. Skiena, Algorithms for square roots of graphs, SIAM J. Discrete Math.,
8 (1995), pp. 99–118.
[90] L. Lovász and M. D. Plummer, Matching Theory, North-Holland, Amsterdam, 1986.
[91] M. Luby, A simple parallel algorithm for the maximal independent set problem, SIAM J.
Comput., 15 (1986), pp. 1036–1053.
[92] B. Manvel, Extremely greedy coloring algorithms, in Graphs and Applications, Proceedingsof the 1st Colorado Symposium on Graph Theory, F. Harary and J. Maybee, eds., John
Wiley, New York, 1985, pp. 257–270.
[93] D. W. Matula, A min-max theorem for graphs with application to graph coloring, SIAM
Rev., 10 (1968), pp. 481–482.
[94] D. W. Matula, G. Marble, and J. Isaacson, Graph coloring algorithms, in Graph Theory
and Computing, R. Read, ed., Academic Press, New York, 1972, pp. 109–122.
[95] D. W. Matula, Y. Shiloach, and R. E. Tarjan, Two Linear-Time Algorithms for FiveColoring a Planar Graph, Tech. Report STAN-CS-80-830, Computer Science Department,
Stanford University, Stanford, CA, 1980.
[96] D. W. Matula, Y. Shiloach, and R. E. Tarjan, Analysis of two linear-time algorithms forﬁve-coloring a planar graph, Congr. Numer., 33 (1981), p. 407 (abstract).
[97] S. T. McCormick, Optimal approximation of sparse Hessians and its equivalence to a graphcoloring problem, Math. Program., 26 (1983), pp. 153–171.
[98] J. Mitchem, Every planar graph has an acyclic 8-coloring, Duke Math. J., 41 (1974), pp. 177–
181.
[99] R. Motwani, Average-case analysis of algorithms for matchings and related problems, J.
ACM, 41 (1994), pp. 1329–1356.
[100] U. Naumann, Cheaper Jacobians by simulated annealing, SIAM J. Optim., 13 (2002), pp. 660–
674.
[101] U. Naumann, Optimal accumulation of Jacobian matrices by elimination methods on the dual
computational graph, Math. Program., 99 (2004), pp. 399–421.
[102] J. Neśetŕil and P. Ossan de Mendez, Colorings and homomorphisms of minor closed
classes, in Discrete and Computational Geometry: The Goodman-Pollack Festschrift,
B. Aronov, S. Basu, J. Pach, and M. Sharir, eds., Springer-Verlag, Berlin, 2003, pp. 651–
664.
[103] G. N. Newsam and J. D. Ramsdell, Estimation of sparse Jacobian matrices, SIAM J.
Algebraic Discrete Methods, 4 (1983), pp. 404–418.
[104] J. J. Nocedal and S. J. Wright, Numerical Optimization, Springer-Verlag, New York, 1999.
[105] W. Nooy, A. Mrvar, and V. Batageli, Exploring Social Network Analysis with Pajek,
Cambridge University Press, Cambridge, UK, 2004.
[106] B. Norris, S. Balay, S. Benson, L. Freitag, P. D. Hovland, L. C. McInnes, and
B. Smith, Parallel components for PDEs and optimization: Some issues and experiences, Parallel Comput., 28 (2002), pp. 1811–1831.
[107] M. J. D. Powell and P. L. Toint, On the estimation of sparse Hessian matrices, SIAM J.
Numer. Anal., 16 (1979), pp. 1060–1074.
[108] P. Raghavan and E. Upfal, Eﬃcient routing in all-optical networks, in Proceedings ofthe 26th Symposium on the Theory of Computing (STOC’94), ACM, New York, 1994,
pp. 134–143.
[109] E. Ramadan, A. Tarafdar, and A. Pothen, A hypergraph model for the yeast protein
complex network, in Proceedings of the IPDPS Workshop on High Performance Computational Biology, CD-ROM, IEEE, 2004.
[110] N. Rostaing, S. Dalmas, and A. Galligo, Automatic diﬀerentiation in Odyssee, Tellus,
45a (1993), pp. 558–568.
[111] Y. Saad, ILUM: A multi-elimination ILU preconditioner for general sparse matrices, SIAM
J. Sci. Comput., 17 (1996), pp. 830–847.
[112] J. Schmidt-Pruzan, E. Shamir, and E. Upfal, Random hypergraph coloring algorithms andthe weak chromatic number, J. Graph Theory, 8 (1985), pp. 347–362.
[113] G. Szekeres and H. S. Wilf, An inequality for the chromatic number of a graph, J. Combin.
Theory, 4 (1968), pp. 1–3.
[114] C. Thomassen, Every planar graph is 5-choosable, J. Combin. Theory Ser. B, 62 (1994),
pp. 180–181.
[115] B. Toft, Coloring, stable sets and perfect graphs, in Handbook of Combinatorics, Vol. I,
R. L. Graham, M. Grötschel, and L. Lovász, eds., MIT Press, Cambridge, UK, 1996,
pp. 233–288.

GRAPH COLORING FOR COMPUTING DERIVATIVES

705[116] M. Trick, Network resources for coloring a graph, http://mat.gsia.cmu.edu/COLOR/
color.html.
[117] V. V. Vazirani, Approximation Algorithms, Springer-Verlag, Berlin, 2001.
[118] A. Walther, ADOL-C: A Package for Automatic Diﬀerentiation of Algorithms Written in
C/C++, http://www.math.tu-dresden.de/wir/project/adolc/.
[119] G. Wegner, Graphs with Given Diameter and a Coloring Problem, Tech. Report, Universityof Dortmund, Germany, 1977.
[120] D. J. A. Welsh and M. B. Powell, An upper bound for the chromatic number of a graphand its application to timetabling problems, Comput. J., 10 (1967), pp. 85–86.
[121] A. Wigderson, Improving the performance guarantee for approximate graph coloring, J.
ACM, 30 (1983), pp. 729–735.
[122] R. A. Wilson, Graphs, Colourings, and the Four-Colour Theorem, Oxford University Press,
Oxford, 2002.
[123] L. A. Wolsey, Integer Programming, Wiley-Interscience, New York, USA, 1998.
[124] M. Yannakakis, Computing the minimum ﬁll-in is NP-complete, SIAM J. Algebraic Discrete
Methods, 2 (1981), pp. 77–79.