INFORMS Journal on Computing

informs

Vol. 21, No. 2, Spring 2009, pp. 209–223
issn 1091-9856  eissn 1526-5528  09  2102  0209®

doi 10.1287/ijoc.1080.0286© 2009 INFORMS

Efﬁcient Computation of Sparse Hessians Using
Coloring and Automatic Differentiation
Assefaw H. Gebremedhin, Arijit Tarafdar

Department of Computer Science and Center for Computational Sciences, Old Dominion University,
Norfolk, Virginia 23529 {agebreme@purdue.edu, tarafdar@cs.odu.edu}

Alex Pothen

Department of Computer Science and Computing Research Institute, Purdue University,
West Lafayette, Indiana 47907, apothen@purdue.edu

Andrea Walther

Institute of Scientiﬁc Computing, Technical University of Dresden, D-01062 Dresden, Germany,
andrea.walther@tu-dresden.deT

he computation of a sparse Hessian matrix H using automatic differentiation (AD) can be made efﬁcient
using the following four-step procedure: (1) Determine the sparsity structure of H, (2) obtain a seed matrix S
that deﬁnes a column partition of H using a specialized coloring on the adjacency graph of H, (3) compute thecompressed Hessian matrix B ≡ HS, and (4) recover the numerical values of the entries of H from B.
The coloring variant used in the second step depends on whether the recovery in the fourth step is direct orindirect: a direct method uses star coloring and an indirect method uses acyclic coloring. In an earlier work, we had
designed and implemented effective heuristic algorithms for these two NP-hard coloring problems. Recently,
we integrated part of the developed software with the AD tool ADOL-C, which has recently acquired a sparsity
detection capability. In this paper, we provide a detailed description and analysis of the recovery algorithms
and experimentally demonstrate the efﬁcacy of the coloring techniques in the overall process of computing the
Hessian of a given function using ADOL-C as an example of an AD tool. We also present new analytical results
on star and acyclic coloring of chordal graphs. The experimental results show that sparsity exploitation via
coloring yields enormous savings in runtime and makes the computation of Hessians of very large size feasible.
The results also show that evaluating a Hessian via an indirect method is often faster than a direct evaluation.
This speedup is achieved without compromising numerical accuracy.
Key words: sparse Hessian computation; acyclic coloring; star coloring; automatic differentiation; combinatorial
scientiﬁc computing
History: Accepted by Karen Aardal, Area Editor for Design and Analysis of Algorithms; received
October 2006; revised August 2007, April 2008; accepted April 2008. Published online in Articles in Advance
September 15, 2008.

1.

Introduction

storage required to compute the Hessian using AD (orestimate it using ﬁnite differences). Matrix compression
has been found to be an effective technique in thiscontext: given an n × n Hessian matrix H of known
sparsity structure, an n × p seed matrix S, with p as
small as possible, is determined; then the numericalvalues in the n × p compressed matrix B ≡ HS are
obtained using AD; and ﬁnally, the nonzero entries ofthe original Hessian H are recovered from the compressed representation B.
One way a matrix H is compressed to form amatrix B is to partition the columns of H into p structurally disjoint groups and let each column of B be thesum of the columns of H that belong to the samegroup. Mathematically, S is an n × p matrix whosej k entry is one if column hj of the matrix H belongs
to group k, and zero otherwise. Because the matrix S
deﬁnes a partitioning of the columns of H, each of its

1.1. Background
Solvers for nonlinear optimization problems suchas IPOPT (Wächter and Biegler 2006) and LOQO
(Vanderbei and Shanno 1999) require second-orderderivatives of the Lagrangian function. Furthermore,exact Hessians are needed in parametric sensitivity analysis, such as in the control of dynamicalsystems in real time (Büskens and Maurer 2001).
A Hessian matrix or Hessian-vector products canbe computed accurately using automatic differentiation(AD), a generic name for a set of chain–rule-based
techniques for evaluating derivatives of a function
given as a computer program (Griewank 2000).
A Hessian matrix that arises in a large-scale application is typically sparse. Sparsity, along with symmetry, can be exploited to reduce the runtime and the
209

210

Gebremedhin et al.: Efﬁcient Computation of Sparse Hessians Using Coloring and AD

entries is either zero or one, and in every row r of S
there exists exactly one column c where the entry src
is equal to one.
The criteria used to ﬁnd a suitable seed matrix S
where the parameter p is minimized—the speciﬁc
partitioning problem—depend on whether the numerical values of the entries of the Hessian H are to be
recovered from the compressed matrix B directly or
indirectly (via substitution). In a direct recovery, nofurther arithmetic is required, whereas in a recovery
via substitution, a set of simple triangular systems ofequations needs to be solved implicitly. The partitioning requirements for a seed matrix suitable for a substitution method are less strict than the requirements
for a seed matrix appropriate for a direct method.
Hence, the former usually results in a smaller p compared with the latter.
The computation of a seed matrix requires only
structural information. Such matrix problems can beformulated conveniently and solved as graph problems. Coleman and Moré (1984) showed that a starcoloring of the adjacency graph of a Hessian models
the partitioning problem that occurs in the computation of the Hessian via a direct method. Coleman
and Cai (1986) showed that the corresponding modelin a substitution-based computation is acyclic coloring. These specialized graph-coloring problems aredescribed in §2.
In a previous work, we had developed effective
heuristic algorithms for the NP-hard star- and acycliccoloring problems, and showed that the algorithms
are superior to previous approaches (Gebremedhin
et al. 2007). In an even earlier work (Gebremedhin
et al. 2005), we had provided a comprehensive reviewof the role of graph coloring in the efﬁcient computation of Jacobians and Hessians. Recently, we
have developed efﬁcient algorithms for recovering thenumerical values of a Hessian from a compressed
representation obtained using a star or an acyclic
coloring. These algorithms take advantage of twocolored structures that the associated coloring algorithms maintain. Our implementations of the coloring
as well as the recovery algorithms, which form
part of the software package ColPack (Gebremedhin
et al. 2009), have been integrated with ADOL-C, an
operator-overloading based AD tool for the differentiation of functions speciﬁed in C/C++ (Griewank
et al. 1996). Walther (2008) recently developed a sparsity detection technique for Hessians and added thefunctionality to ADOL-C.
1.2. Contributions
Our contributions in this paper are analytical as well
as experimental. In terms of theory, in §3 we provide
a detailed description and analysis of the Hessianrecovery algorithms mentioned earlier. Our analysis

INFORMS Journal on Computing 21(2), pp. 209–223, © 2009 INFORMS

includes time complexity as well as numerical stability: we show that our recovery algorithms are linear
in the size of the problem and that indirect recovery
using two-colored trees is nearly as stable as direct
recovery. We also present, in §4, new results on star
and acyclic coloring of chordal graphs, a class towhich adjacency graphs of banded matrices belong.
Experimentally, we demonstrate in §5 the advantage that coloring techniques offer in sparse Hessiancomputation using AD. We use two test cases:
a real-world power transmission problem and a synthetic unconstrained quadratic optimization problem. In both cases, we compute Hessians of various
dimensions and sparsity structures, including banded
and random structures, using ADOL-C. The results
obtained show that sparsity exploitation via star andacyclic coloring enables one to affordably compute
Hessians of dimensions that could not have been
computed otherwise. For sizes where dense Hessiancomputation is at least possible, the saving in runtime obtained by exploiting sparsity via coloring is
dramatic. Furthermore, an indirect method that uses
acyclic coloring quite often is found to be faster than adirect method that relies on star coloring, considering
the overall process. This speedup is achieved without compromising numerical accuracy. The experimental results also show that both the star- and theacyclic-coloring heuristic algorithm we used ﬁnd optimal solutions for banded matrices.
1.3. Related WorkIn sparse Hessian computation, an approach that canbe used either orthogonal to or in combination withcompression is the use of elimination techniques on thecomputational graph of the Hessian. This approach
was ﬁrst considered by Dixon (1991) and is the subjectof current research in the AD community. Partial separability is another related approach used in sparse
Hessian computation. Gay (1996) has considered suchan approach and its implementation using AMPL.

2.

Graph Coloring

In the current work, we have used the star- andacyclic-coloring algorithms we had developed earlier(Gebremedhin et al. 2007). After reviewing a few preliminary concepts, we brieﬂy discuss these algorithms
in this section.
2.1. Preliminaries
Two distinct vertices in a graph are distance-k neighbors if a shortest path connecting them consists ofat most k edges. The degree-k of a vertex v, denotedby dk v , is the number of distinct paths of length
at most k edges that start at the vertex v; the average
degree-k is denoted by d¯k v . A distance-k coloring ofa graph G = V  E is a mapping  V → 1 2     pGebremedhin et al.: Efﬁcient Computation of Sparse Hessians Using Coloring and AD
INFORMS Journal on Computing 21(2), pp. 209–223, © 2009 INFORMS

such that whenever vertices v and w are distance-k
neighbors, v = w . The kth power of a graphG = V  E is the graph Gk = V  F , where v w ∈ F
whenever vertices v and w are distance-k neighbors
in G. A mapping  is a distance-k coloring of agraph G if and only if  is a distance-1 coloring of
the graph Gk .
A star coloring of a graph is a distance-1 coloring
where, in addition, every path in the graph on four
vertices (P4 ) is required to use at least three colors.
An acyclic coloring of a graph is a distance-1 coloring with the further restriction that every cycle in
the graph uses at least three colors. The names star
and acyclic coloring are due to the structure of twocolored induced subgraphs. In a star-colored graph, asubgraph induced by any two color classes—sets ofvertices having the same color—is a collection of stars.
A star is a complete bipartite graph in which one ofthe vertex sets consists of a single vertex, called thehub. The other vertices of the star are spokes. An edge
is a degenerate case of a star, in which one vertex is
the hub and the other, the spoke, assigned arbitrarily.
Similarly, in an acyclically colored graph, a subgraphinduced by any two color classes is a collection oftrees, and thus acyclic.
Each one of the distance-k, star-, and acycliccoloring problems, whose objective is to appropriately color a graph using the fewest possible colors, is known to be NP-hard (Coleman and Moré
1984, Coleman and Cai 1986, Lin and Skiena 1995).
The minimum number in each case is referred toas a specialized chromatic number of the graph. The
distance-k chromatic number k G , the star chromatic number s G , and the acyclic chromatic number a G of a general graph G are all known tobe hard to approximate (Gebremedhin et al. 2007).
The following observation on the relationship among
these chromatic numbers follows directly from thedeﬁnitions of the various colorings and the fact thata distance-1 coloring of a graph G is equivalent to a
distance-2 coloring of the square graph G2 .
Observation 2.1. For every graph G, 1 G ≤a G ≤ s G ≤ 2 G = 1 G2 .
2.2. Star- and Acyclic-Coloring Algorithms
The star- and acyclic-coloring heuristic algorithms wedeveloped in Gebremedhin et al. (2007) are founded
on one common key idea: maintain and efﬁciently use
the structure of two-colored induced subgraphs. As
stated earlier, the respective subgraphs here are a collection of stars or trees. Both algorithms are greedy
in the sense that a partial coloring is progressively
extended by processing one vertex at a time. In each
step, a vertex v of an input graph G = V  E is
assigned the smallest allowable color with respect tothe current partial coloring. This is achieved by ﬁrst

211

determining the set of forbidden colors to the vertex v. Here, identifying and forbidding colors used by
the vertices adjacent to the vertex v is straightforward.
Maintaining the collection of two-colored structures
turns out to be crucial for identifying and forbidding colors that could lead to two-colored paths onfour vertices in the case of star coloring, and to twocolored cycles in the case of acyclic coloring.
For star coloring, there are exactly two cases thatneed to be considered to avoid a two-colored P4
involving the vertex v.
1. If the vertex v has two adjacent vertices w and xof the same color, then the color of every vertex adjacent to w and the color of every vertex adjacent to x
should be marked as a forbidden color to the vertex v.
2. If the vertex v has exactly one adjacent vertex w
of color c, and w belongs to a star S consisting ofat least two edges but is not the hub of S, then thecolor used by the hub of S should be marked as aforbidden color to the vertex v.
Likewise, for acyclic coloring, there is one case thatneeds to be addressed to avoid a two-colored cycle
involving the vertex v:
1. If the vertex v is adjacent to two vertices in a
two-colored tree T , and these vertices have the samecolor, then the other color used in the tree T should
be marked as a forbidden color to the vertex v.
Once the color for the vertex v is determined, everyedge incident on v that leads to an already-coloredvertex is placed in the appropriate two-colored structure. Note that because a star (respectively, an acyclic)coloring is also a distance-1 coloring, there is a unique
two-colored star (respectively, two-colored tree) towhich such an edge belongs. In other words, a star
(respectively, an acyclic) coloring partitions the edge
set of a graph into two-colored stars (respectively,
two-colored trees). In the acyclic-coloring algorithm
being described, the placement of an edge in a twocolored tree may result in the merging of two previously disconnected trees in a two-colored forest. Sucha situation does not arise in the star-coloring case.
In Gebremedhin et al. (2007) it was shown thatthe acyclic-coloring algorithm sketched above can be
implemented efﬁciently using the disjoint-set data
structure to maintain the collection of two-colored
trees. A simpler data structure was shown to be sufﬁcient for a similar purpose in the star-coloringalgorithm. The complexity of the star- and theacyclic-coloring algorithm was shown to be OV d¯2
and OV d¯2 · , respectively, where is the inverseof Ackermann’s function, a function that grows very
slowly with V .

3.

Sparse Hessian Computation

We now relate the graph-coloring variants discussed
in the previous section with partitioning strategies for
Hessian compression for direct and indirect recovery.

Gebremedhin et al.: Efﬁcient Computation of Sparse Hessians Using Coloring and AD

212

INFORMS Journal on Computing 21(2), pp. 209–223, © 2009 INFORMS

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



h1

h2
1

h3
2

1

1

4
h6

5

1

h5

3
h7

h11

 h21 + h23 + h25



h33



h43



h55


B1 = HS1 =  h63 + h65 + h69



h71


 h85 + h89



h99


h109

h4

h8

3

2

h9

h12

h17

0

h22

0

0

h32

h34

h36

h410

h44

0

h52

0

h56

0

0

h66

0

h77

0

0

h87

0

h910

0

h96

h1010

h104

0

h10

0




0 


0 


0 


h58 



0 


h78 


h88 


h98 


0

Figure 1

Top—A Symmetrically Orthogonal Partition of the Columnsof a Hessian and Its Representation as a Star Coloring of the
Adjacency Graph; Bottom—The Compressed Matrix Obtainedby Adding Together Columns that Belong to the Same Groupin the Partitioning
Note. The illustration uses ﬁve colors—red, blue, yellow, green, and navy
blue—that are also represented by the integers 1, 2, 3, 4, and 5, respectively,
as shown at the bottom edge of the matrix as well as adjacent to the verticesof the graph.

3.1.

Direct Method

3.1.1. Matrix Partitioning. For a direct recovery,
the columns of a Hessian matrix H need to be partitioned in such a way that for every nonzero entry
hij , either hij or its symmetric counterpart hji appears
as a sole entry in the compressed matrix. A symmetrically orthogonal partition of H precisely captures
this requirement. A partition of the columns of a Hessian H is symmetrically orthogonal if for every nonzero
element, hij , either (1) the group containing column
hj has no other column with a nonzero in row i, or(2) the group containing column hi has no other column with a nonzero in row j.
The left part of the upper row of Figure 1 illustrates
a symmetrically orthogonal partition of a Hessian
matrix H with a speciﬁc sparsity pattern. A nonzero
entry of H is denoted by the symbol “X” and a zero
entry is left blank. The 10 columns of H are partitioned into ﬁve groups, and columns that belong
to the same group are painted with the same color.
The ﬁve groups in the partition—represented by thecolors red, blue, yellow, green, and navy blue—are
also identiﬁed by the integers 1 through 5, respec-tively, as shown at the bottom edge of the matrix.
(See the Online Supplement to this article, available
at http://joc.pubs.informs.org/ecompanion.html, fora full-color version of this illustration, as well as theone for Figure 4.) The partition in the illustration in
Figure 1 deﬁnes a 10×5 seed matrix, which we denote
by S1 . For example, column 1 of S1 has 1s in rows 1,
3, 5, and 9, corresponding to the columns of H thatbelong to group 1 (color red), and 0s in all other rows.
The lower row of Figure 1 shows the resultant compressed matrix B1 = HS1 . The reader can easily verify
that every nonzero entry of the matrix H (or its symmetric counterpart) can be read off directly from someentry of the matrix B1 .
3.1.2. Coloring Model. Let H be a Hessian, eachof whose diagonal elements is nonzero. The adjacency
graph GH of the matrix H is an undirected graphwhose vertex set consists of the columns of H andwhose edge set consists of pairs hi  hj whenever thematrix entry hij , i = j, is nonzero. In such a graph,entries hij and hji are represented by the single edgehi  hj , and there are no explicit edges representing
the diagonal entries of H . Coleman and Moré (1984)
established that the problem of ﬁnding a symmetrically orthogonal partition of a Hessian having the
fewest groups is equivalent to the star-coloring problem on its adjacency graph. The right part of theupper row of Figure 1 illustrates this equivalence.
3.1.3. Recovery Routines. Figure 2 outlines a simple routine, called directRecover1, for recovering thenumerical values of the nonzero entries of a Hessian
H from its compressed representation B obtained via
a star coloring of GH . The routine achieves this byconsidering the structure of H one row at a time.
The nonzero entries in a speciﬁc row are in turnconsidered one element at a time. For each nonzero
entry hij in row i, the if-test in the inner for-loop
checks whether there exists a column index j where
entries hij and hij “collide”; i.e., columns hj and hj
belong to the same group. Depending on the outcome of the test, either the entry hij or the entry hji
Input: The adjacency graph GH of a Hessian H of order n; avertex-indexed integer array color specifying a star coloring ofGH ; a compressed matrix B representing H .
Output: Numerical values in H .
for i ← 1 2     n
for each j where hij = 0
if ∃j = j where hij = 0 and color"hj # = color"hj #H "j i# ← B"j color"hi ##;else
H "i j# ← B"i color"hj ##;end-ifend-forend-forFigure 2

directRecover1—A Routine for Recovering a Hessian froma Star-Coloring-Based Compressed Representation

Gebremedhin et al.: Efﬁcient Computation of Sparse Hessians Using Coloring and AD

213

INFORMS Journal on Computing 21(2), pp. 209–223, © 2009 INFORMS

Input: The adjacency graph GH of a Hessian H of order n; avertex-indexed integer array color specifying a star coloringof GH ; a set  of two-colored stars; a compressed matrix B
representing H .
Output: Numerical values in H .
for i ← 1 2     n
H "i i# ← B"i color"hi ##;end-forfor each two-colored star S ∈ 
Let hj be the hub vertex in S;
for each spoke vertex hi ∈ S
H "i j# ← B"i color"hj ##;end-forend-forFigure 3

directRecover2—A Routine Using Two-Colored Stars for
Recovering a Hessian from a Star-Coloring-BasedCompressed Representationis read from an appropriate location in the matrix B.
Clearly, the if-test can be performed within Od1 hitime, where hi is the vertex being considered in thecurrent iteration of the outer for-loop, and d1 hi is
its degree-1 in the adjacency graph GH . Hence, thecomplexity of directRecover1 is OE · d¯1 , where d¯1is the average vertex degree in GH .
The recovery of Hessian entries in a direct method
theoretically could be done more efﬁciently by using
the set of two-colored stars deﬁned by a star coloring of the adjacency graph. Speciﬁcally, directRecover2, the routine speciﬁed in Figure 3, shows thatthe recovery of the entries of the matrix H fromthe matrix B can be done in OE -time when thetwo-colored structures are readily available. As canbe seen in the ﬁrst for-loop, because adjacent vertices in a star coloring receive different colors, eachdiagonal element hii can be retrieved simply from
B"i colorhi #. The second (outer) for-loop shows thateach off-diagonal element hij can be obtained by consulting the unique two-colored star to which the edgehi  hj belongs. The reader is encouraged to see theroutine in Figure 3 in conjunction with the illustration in Figure 1. For example, one can see that all of
the edges (off-diagonal nonzeros) that belong to the
red-blue (color 1-color 2) star induced by the vertices
h1  h2  h3  h5  can be obtained from the group that
corresponds to the color of the hub vertex h2 (i.e., column 2 of B1 ). Note also that an edge such as h1  h7
that belongs to a single-edge star can be obtained
from either one of its endpoints.
3.2.

Substitution Method

3.2.1. Matrix Partitioning. If the requirement thatthe entries of a Hessian be recovered directly froma compressed representation is relaxed, then thecompression can be done much more compactly. One
possibility here is to use what is called a substitutablepartition. A partition of the columns of a symmetricmatrix H is substitutable if there exists an ordering on
the elements of H such that for every nonzero element
hij , either (1) column hj is in a group where all thenonzeros in row i, from other columns in the samegroup, are ordered before hij ; or (2) column hi is in agroup where all the nonzeros in row j, from othercolumns in the same group, are ordered before hij .
In the deﬁnition just stated, a nonzero entry hij of asymmetric matrix H is identiﬁed with the entry hji .
A nonzero entry hij is said to be ordered before anonzero entry hi j if hij is evaluated before hi j .
3.2.2. Coloring Model. Fortunately, the rather
clumsy notion of substitutable partition has a simple
graph-coloring formulation: Coleman and Cai (1986)
proved that an acyclic coloring of the adjacency graphof a Hessian induces a substitutable partition of its
columns. Thus, the problem of ﬁnding a substitutablepartition with the fewest groups reduces to the problem of ﬁnding an acyclic coloring with the fewest
colors.
We use the illustration in Figure 4 to show how anacyclic coloring can be used to obtain a substitutablepartition. The matrix H shown in Figure 4 has thesame sparsity structure as the matrix in Figure 1. The
1 2 3 4 5 6 7 8 9 10

XX
X
XXX X
XXX X
XX
X
X
XX X
H
X XX
X
X
XX
X X XX
X X XX
X
XX
1 2 1 2 1 3 2 3 2 1

h1

h2
1

h3
2

1

1

3

2

3

2

h8

h9

h11

 h21 + h23 + h25



h33


 h +h

43
410


h55


B2 = HS2 =  h63 + h65



h71



h85



h
910


h1010h12 + h17

0

Figure 4

2

h6

h5

h7



h4

h22
h32 + h34
h44
h52
h69
h77
h87 + h89h99h104 + h109

1
h10







h36 



0


h56 + h58 



h66 


h78 


h88 


h96 + h98 00

Top—A Substitutable Partition of the Columns of a Hessianand Its Representation as an Acyclic Coloring of the
Adjacency Graph; Bottom—The Compressed Matrix Obtainedby Adding Together Columns that Belong to the Same Groupin the Partitioning
Note. The illustration uses three colors—red, blue, and yellow—that are also
represented by the integers 1, 2, and 3, respectively, as shown at the bottomedge of the matrix as well as adjacent to the vertices of the graph.

214

Gebremedhin et al.: Efﬁcient Computation of Sparse Hessians Using Coloring and AD

coloring of the adjacency graph GH depicted in Figure 4 is clearly acyclic. In the illustration, the columnsof the matrix H have been painted in accordance withthe shown acyclic coloring of GH . The lower row
shows the compressed matrix obtained using the seedmatrix S2 deﬁned by the depicted acyclic coloring.
3.2.3. Recovery Routine. We proceed to show
how the acyclic coloring in Figure 4 will be used inrecovering the entries of the matrix H from the compressed matrix B2 = HS2 .
First, observe that every diagonal element hii isdirectly recoverable from the compressed matrix B2 ,because hii appears “alone” in row i and column kof the matrix B2 , where k corresponds to the groupto which column hi belongs. This is a direct consequence of the fact that adjacent vertices in an acyclically colored graph, as in a star-colored graph, receive
different colors.
Next, consider the task of determining the offdiagonal nonzero entries hij , i = j. Recall that eachsuch matrix entry and its symmetric counterpart correspond to an edge in the adjacency graph, and an
acyclic coloring partitions the edges of the graph into
two-colored trees. We use each of these two-colored
trees, separately, to deﬁne an order in which edges canbe solved for.
In a two-colored tree, every edge incident on a leaf
vertex can be determined directly because it is the
only nonzero in a row of the group of columns towhich its parent vertex belongs. As an example, consider the red-blue (color 1-color 2) tree in Figure 4
induced by the vertices h1  h2  h3  h4  h5  h7  h9  h10 .
In this tree, the vertex h7 is a leaf, and the edge h7  h1
can be immediately read from row 7 of the ﬁrst column of B2 , the group to which column h1 belongs.
Similarly, the vertex h9 is a leaf, and the edge h9  h10
can be read from row 9 of the ﬁrst column of B2 , again
the group to which column h10 belongs. Likewise,edge h5  h2 can be directly obtained from B2 "5 2#.
Once the edges incident on leaf vertices have been
determined, they can be deleted from the tree to create new leaves. The process can then be repeated tosolve for edges incident on the new leaf vertices byusing values computed for the leaf edges from earlier
steps. The process terminates when the tree becomes
empty, i.e., when all of the edges have been evaluated. In general, there are alternative ways in which
an edge can be solved for, so the evaluation process
is not unique.
Returning to our illustration, in the red-blue tree,
once the edges h7  h1 , h9  h10 , and h5  h2 havebeen evaluated and deleted, the path h1 –h2 –h3 –h4 –h10
remains. In this path, the edges h1  h2 and h10  h4
are incident on leaf vertices; the edge h1  h2 can beevaluated using B2 "1 2# and the previously computed
value for the edge h7  h1 , and the edge h10  h4 can

INFORMS Journal on Computing 21(2), pp. 209–223, © 2009 INFORMS

Input: The adjacency graph GH of a Hessian H of order n; avertex-indexed integer array color specifying an acyclic coloringof GH ; a set  of two-colored trees; a compressed matrix B
representing H .
Output: Numerical values in H .
for i ← 1 2     n
H "i i# ← B"i color"hi ##;end-forfor each two-colored tree T ∈ 
for each vertex hj ∈ TstoredValues"hj # ← 0;
end-for
while the tree T is not empty
Pick a leaf vertex hi ∈ T ;Let hj be the parent of hi in T ;H "i j# ← B"i color"hj ## − storedValues"hi #;storedValues"hj # ← storedValues"hj # + H "i j#;
Delete vertex hi (along with edge hi  hj ) from T ;end-while
end-forFigure 5indirectRecover—A Routine Using Two-Colored Trees for
Recovering a Hessian from an Acyclic-Coloring-BasedCompressed Representation

be evaluated using B2 "10 2# and the previously computed value for the edge h9  h10 . After this, the edges
h1  h2 and h4  h10 can be deleted, leaving the pathh2 –h3 –h4 from which the edges h2  h3 and h3  h4
can be evaluated.
The red-blue tree in Figure 4 enabled the determination of 7 of the 13 distinct off-diagonal nonzero entries.
The remaining six nonzeros are determined using theother two trees, the red-yellow (color 1-color 3) tree
and the blue-yellow (color 2-color 3) tree.
We summarize the process we have been describing thus far in Figure 5, where we outline the routine
indirectRecover for evaluating the nonzeros of a
Hessian from a compressed representation induced byan acyclic coloring. Note the resemblance between theroutines directRecover2 and indirectRecover. The
ﬁrst for-loops in each case correspond to the determination of diagonal nonzeros, and the second forloops to the recovery of off-diagonal nonzeros fromtwo-colored stars (respectively, trees). In indirectRecover, the variable storedValues is used to store edge
values that will be “substituted” in the determinationof edges in later steps. Speciﬁcally, let hi  hj be a pairof child and parent vertices, respectively, in a twocolored tree T , and let T hi denote the subtree of T
rooted at the vertex hi that would remain if the edgehi  hj were to be removed from T . Then it is easy to
see that the quantity stored in the variable storedValues at the index hi isstoredValues"hi # =H "r s#(1)
hr  hs ∈ET hiwhere ET hi denotes the set of edges in the tree T hi ,and the entry H "i j# of the Hessian is computed usingGebremedhin et al.: Efﬁcient Computation of Sparse Hessians Using Coloring and AD

215

INFORMS Journal on Computing 21(2), pp. 209–223, © 2009 INFORMS

the equation
H"i j# = B"i color"hj ## − storedValues"hi #

(2)

Using appropriate data structures, the computational
work associated with each two-colored tree in indirectRecover can be performed in time proportional
to the number of edges in the tree. Thus, the overall complexity of indirectRecover on an adjacency
graph GH = V  E is OE .
3.3. Numerical Stability
How does the routine indirectRecover compare with the routines directRecover1 or directRecover2 in terms of numerical stability? The answer
to this question turns out to be quite positive:
indirectRecover for practical purposes is nearly as
stable as directRecover. Our analytical justiﬁcation
for this claim is a natural extension of the works
of Powell and Toint (1979) and Coleman and Cai(1986), who analyzed the error bounds associatedwith (specialized) substitution methods in the context
of Hessian estimation using ﬁnite differences. In ourcontext, because the compressed Hessian B is computed analytically using automatic differentiation (andthus exactly within machine precision), the associatednumerical stability analysis is fundamentally different. Issues such as truncation error and choice of step
length do not arise in our context.
The only arithmetic operations involved in indirectRecover are subtraction and addition; theabsence of division is highly favorable for numericalstability. Furthermore, the determination of nonzeros(edges) in one two-colored tree is entirely independent of the determination of edges in another twocolored tree, and therefore there is less opportunity
for error magniﬁcation. As we shall show shortly, theerror accumulation within a two-colored tree is inturn very limited.
Let hi  hj be an edge in the input graph GH =V  E to indirectRecover, and let T = VT  ET
be the two-colored tree to which the edge hi  hj
belongs. Furthermore, as done earlier, let T hi =VT hi  ET hi denote the subtree of T rooted at the vertex hi that is obtained by removing the edge hi  hjfrom T . Let n, nT , and nT hi denote the number ofvertices in GH , T , and T hi , respectively. Our goal
is to prove a bound on the accuracy of the Hessian
entry H"i j# computed using Equation (2). First, following Powell and Toint (1979), we deﬁne the error
matrix , the pointwise difference between the computed Hessian H and the analytic Hessian $ 2 f , as2

"i j# = H "i j# − $ f x

ij

(3)

for all pairs i and j such that hi  hj is an edge inGH . In Theorem 3.1 we will show that "i j#, themagnitude of the error associated with the computation of H "i j# by indirectRecover, is bounded by the
product of nT hi , the number of vertices in the subtree
T hi of T , and a constant independent of T .
Theorem 3.1. The numerical value computed by indirectRecover for each edge hi  hj in the input graphGH is such that "i j# ≤ nT hi · &, where & is a positive
constant.
Proof. Because the compressed Hessian B is computed in ﬂoating point arithmetic, it inevitably contains rounding errors. LetB"icolor"hj ## = B"icolor"hj ##+'"icolor"hj ##(4)denote the computed matrix taking such errors into
account. Thus, the actual value H "i j# evaluated using
indirectRecover is color"hj ## − storedValues"hi #
H "i j# = B"i

(5)

Analogous to the error matrix  associated with H , let
the error matrix ( be the pointwise difference between
the computed values contained in B and the corresponding values in the analytic Hessian $ 2 f ; i.e., let("icolor"hj ##= B"icolor"h
j ##−



$ 2 f xhr hs ∈ET hirs −$

2

f x

ij 

(6)

Using Equations (1) and (5), Equation (6) can be written as("i color"hj ## = H "i j# − $ 2 f x ij+
H "r s# −  $ 2 f x
hr  hs ∈ET hi= "i j# +rs

"r s#hr  hs ∈ET hi

It then follows that
"i j# = ("i color"hj ## −"r s#hr  hs ∈ET hi

Applying the same decomposition to each "r s# forhr  hs ∈ ET hi , one obtains
"i j# = ("i color"hj ## −("r color"hs ##
hr  hs ∈ET hiTaking the absolute values of scalar quantities andnoting that the tree T hi has nT hi − 1 edges,
"i j# ≤ ("i color"hj ##+
("r color"hs ##

(7)hr  hs ∈ET hi≤ ("i color"hj ## + nT hi − 1·maxhr  hs ∈ET hi

("r color"hs ##(8)

Gebremedhin et al.: Efﬁcient Computation of Sparse Hessians Using Coloring and AD

216

INFORMS Journal on Computing 21(2), pp. 209–223, © 2009 INFORMS

If we assume that the second derivatives of f arebounded, then there exists a positive constant M suchthat the maximum of the two terms in Equation (8),and indeed of similar terms in the entire graph GH ,can be bounded as follows:& = max ("r  color"hs ##1≤r  s ≤n

≤ M + max '"r  color"hs ##
1≤r  s ≤nThus, Equation (8) reduces to "i j# ≤ nT hi · &,which is what we wanted to show. 
Suppose the acyclic coloring used in the context
of indirectRecover is actually a star coloring. Then,clearly, for each edge hi  hj , the two-colored tree towhich the edge hi  hj belongs is a star; i.e., T hi issimply the vertex hi . In such a case, in agreement with
our expectations, Theorem 3.1 suggests that the error
associated with the evaluation of hi  hj is bounded
simply by &.

4.

Coloring Chordal Graphs

The test suite in our experiments includes Hessian
matrices with banded nonzero structures, whose adjacency graphs are band graphs (deﬁned later in thissection). Here, we present analytical results on the
distance-k, star, and acyclic chromatic numbers of thelarger class of chordal graphs, an important class ofgraphs with a wide range of applications (Brandstädt
et al. 1999).
Let A be a symmetric matrix of order n with
nonzero diagonal elements. The lower bandwidth of A
is deﬁned as +l A = maxi − j i > j aij = 0, and the
bandwidth of A is the quantity +A = 2+l A + 1. The
matrix A is banded if it is completely dense within
the band; i.e., if for any pair of indices 1 ≤ i, j ≤ n,
i − j ≤ +A ⇔ aij = 0.
The bandwidth of a matrix has a twin concept in
the adjacency graph.
Let G = V  E be a graph on n vertices and let .
be an ordering v1  v2      vn of the vertices. The bandwidth of the ordering . in G is +. G = maxi − j 
vi  vj ∈ E, and the bandwidth of the graph G is +G =min+. G  . is an ordering of V . A graph G is aband graph if there exists an ordering v1  v2      vn of
its vertices such that for any pair of indices i j drawn
from the set 1 2     n, i − j ≤ +G ⇔ vi  vj ∈ E;
the order v1  v2      vn is referred to as the natural ordering of the band graph G. For a generalgraph, ﬁnding a vertex ordering with the minimum
bandwidth—computing the quantity +G —is an NPcomplete problem (Papadimitriou 1976).
The bandwidth of a symmetric matrix A and that of
the adjacency graph GA are related in the following
way: +GA ≤ +l A ≡ +A − 1 /2. In this relationship, equality holds when A is banded, in which caseGA is a band graph, and the natural ordering of
the vertices of GA corresponds to the given ordering of the columns of A. When A is not banded, there
exists a permutation matrix P such that +GA =+l PAP T .
A graph G is chordal if every cycle in G of length
at least four has a chord—an edge that connects two
nonconsecutive vertices in the cycle. Clearly, a bandgraph is chordal. In fact, it is a highly structured,almost regular, chordal graph: the degree dvi of the
ith vertex in the natural ordering of the vertices can be
expressed as dvi = dvn−i+1 = +G + i − 1 for 1 ≤ i ≤+G , and dvi = 2+G for +G + 1 ≤ i ≤ n − +G .
A graph does not need to be this regular to be chordal.
For example, the adjacency graph of a symmetricmatrix in which rows are allowed to have variable
number of nonzeros, but the nonzeros in every row
are required to be consecutive, is chordal but notband.
In what follows we present results, which we
believe to be new, concerning the relationships among
the distance-k chromatic number k G , the star chromatic number s G , the acyclic chromatic number a G , the clique number 1G —the size of thelargest clique in G—and the bandwidth +G of achordal (not necessarily band) graph G. We begin
with a simple observation that is true of any (not necessarily chordal) graph.
Lemma 4.1. For every graph G = V  E , 1G ≤+G + 1, and 1G2 ≤ min2+G + 1 V .
Proof. Let . be an ordering of the vertices of G
such that +. G = +G . Suppose there exists a clique
Q in G of size greater than +G + 1. Then it meansthat there exists some pair of vertices v and w in
the clique Q such that .v − .w  > +G , a contradiction. Hence, the clique number of G cannot
exceed +G + 1. The result for the square graph canbe shown in an analogous fashion. 
Lemma 4.2. If a mapping  is a distance-1 coloring fora chordal graph G, then  is also an acyclic coloring for G.
Proof. Let  be a distance-1 coloring for a chordal
graph G. Consider any cycle C in G, and let l be its
length. If l = 3, then  clearly uses three colors. Ifl ≥ 4, then C contains a chord, and therefore  uses
at least three colors. Hence,  is an acyclic coloringfor G. 
Theorem 4.3. For every chordal graph G, a G =
1 G = 1G ≤ +G + 1. In the last relationship, equality
holds when G is a band graph.
Proof. The equality a G = 1 G follows from
Lemma 4.2. Because a chordal graph is perfect, 1 G =1G by the perfect graph theorem (Lovász 1972). The
last inequality was proven in Lemma 4.1 for anyGebremedhin et al.: Efﬁcient Computation of Sparse Hessians Using Coloring and AD
INFORMS Journal on Computing 21(2), pp. 209–223, © 2009 INFORMS(including chordal) graph. When G is a band graph,clearly, 1G = +G + 1. 
There exist chordal graphs where the last inequality
in Theorem 4.3 is strict. For example, a star graph G
on n vertices has 1G = 2, but +G = n/2.
Theorem 4.4. For every chordal graph G = V  E ,s G ≤ 2 G = 1G2 ≤ min2+G + 1 V 
In both the ﬁrst and the third relationship, equality holds
when G is a band graph.
Proof. The ﬁrst inequality follows from Observation 2.1, and the third from Lemma 4.1. The squaregraph of a chordal graph is chordal and hence perfect.
Therefore, 2 G = 1 G2 = 1G2 . Turning to the special case of band graphs, Coleman and Moré (1984)
have shown that for a band graph G with at least
3+G + 1 vertices, s G = 2 G . For a band graph G,
1G2 = min2+G + 1 V . 
There exist chordal graphs where the ﬁrst inequality in Theorem 4.4 is strict. An example, once again,
is a star graph G on n vertices, which has s G = 2
but 2 G = n.
If symmetry were to be ignored, a structurally orthogonal partition of the columns of a Hessian—a partitionin which no two columns in a group have nonzeros
at a common row index—could be used to compress
a Hessian in a direct method. As McCormick (1983)
ﬁrst showed, a structurally orthogonal partition of a
Hessian can be modeled by a distance-2 coloring of
its adjacency graph. In light of these facts, the results G = 2 G = 2+G + 1 for a band graph G given in
Theorem 4.4 is negative: it shows that exploiting symmetry in a direct computation of a banded Hessian
matrix (star coloring) does not lead to fewer colors
in comparison with a direct computation that ignores
symmetry (distance-2 coloring). The result 1 G =a G = +G + 1 in Theorem 4.3, on the contrary,shows that symmetry exploitation in a banded matrix
is worthwhile in a substitution method.
We conclude this section with some remarks on theperformance of greedy coloring algorithms on chordal
graphs. Recall that a greedy coloring algorithm processes vertices in some order, each time assigninga vertex the smallest allowable color subject to theconditions of the speciﬁc coloring problem. Thereexist several “degree”-based ordering techniques—
including largest-degree-ﬁrst, smallest-degree-last andincidence-degree—that have proven to be quite effective (but still suboptimal) for distance-k coloring ofgeneral graphs (Gebremedhin et al. 2005).
For chordal graphs, better ordering techniques
exist. Given a graph G = V  E , an orderingv1  v2      vn of the vertices in V is a perfect elimination ordering (peo) of G if for all i ∈ 1 2     n, the217vertex vi is such that its neighbors in the subgraphinduced by the set vi      vn  form a clique. It is wellknown that a graph G is chordal if and only if it has apeo. It is also known that a greedy distance-1 coloringalgorithm that uses the reverse of a peo of G gives anoptimal solution, i.e., computes a coloring with 1 G
colors (Brandstädt et al. 1999). For the special case of
a band graph, the natural ordering of the vertices, aswell as its reverse, is a peo. Thus, a greedy distance-1
coloring algorithm that uses the natural ordering of
the vertices would give an optimal solution. However,as Lemma 4.5 and its corollary will imply, an optimalcoloring for a band graph can be obtained without
actually executing the greedy algorithm.
Lemma 4.5. Let G = V  E be any graph and let
v1  v2      vn be an ordering in which the bandwidth of G
is attained. Then the mappings 1 vi = i mod +G + 1
and 2 vi = i mod 2+G + 1 deﬁne a distance-1 coloring and a distance-2 coloring of G, respectively. If G is aband graph, then both of these colorings are optimal.
Corollary 4.6. For every graph G = V  E , 1 G ≤+G + 1, and 2 G ≤ min2+G + 1 V . In both relationships, equality holds when G is a band graph.
The optimality of 1 and 2 in Lemma 4.5 in thecase of band graphs, and the implied equalities inCorollary 4.6, follow from Theorems 4.3 and 4.4.

5.

Numerical Experiments

In this section, we present experimental results concerning the following four steps involved in the
efﬁcient computation of the Hessian of a givenfunction f .
Step S0: Detect the sparsity pattern of the Hessian
Step S1: Obtain a seed matrix S using an appropriate graph coloring
Step S2: Compute the Hessian-seed matrix productB ≡ HS
Step S3: Recover the nonzero entries of H from B
We use two different optimization problems as test
cases: an electric power ﬂow problem, representing a real-world application; and an unconstrained
quadratic optimization problem, a synthetic case chosen for a detailed performance analysis. The underlying test function f in both test cases is speciﬁed in the
programming language C. The coloring and recovery
codes (steps S1 and S3) are written in C++ as partof the software package ColPack (Gebremedhin et al.
2009) and are incorporated into ADOL-C. For step S3,
in the direct case, the routine directRecover1 isused. For the step S2, the second-order adjoint modein the latest version of ADOL-C, which is signiﬁcantly faster than previous versions, is used. Whenneeded, the sparsity structure of the Hessian (step S0)

Gebremedhin et al.: Efﬁcient Computation of Sparse Hessians Using Coloring and AD

218

INFORMS Journal on Computing 21(2), pp. 209–223, © 2009 INFORMS

Table 1

Absolute Runtimes in Seconds for the Evaluation of the Function f and the Steps S0, S1, S2, and S3 for Test Hessiansin the Optimal Power Flow Problem
DirectIndirectnEvalf 

S0

S1

S2

S3

Tot

S1

S2

S3

Tot

Dense

72
3,760
4,472
9,932
22,5400.000008
0.000587
0.000717
0.001933
0.002118

0001
0247
0372
6212
226620.00060.02520.03290.09290.2170

0007
0805
0899
2839
201620.00003
0.002480.00307
0.007190.01741

0008
0833
0935
2939
203960.00040.04640.06010.22201.0402

0005
0646
0561
1159
125310.00040.06170.07910.30491.2981

0005
0755
0700
1686
14869008
2013225799131543∗∗∗

Notes. The last column shows runtime for the computation of a Hessian without exploiting sparsity. The asterisks indicate that space
could not be allocated for the full Hessian.

is determined using the recently added functionality in ADOL-C (Walther 2008). The experiments on
the power ﬂow problem are performed on a Linuxsystem with an Intel Xenon 1.5 GHz processor and
1 GB RAM, and those on the synthetic problem are
performed on a Fedora Linux system with an AMD
Athlon XP 1.666 GHz processor and 512 MB RAM. In
both cases, the gcc 4.1.1 compiler with -02 optimization is used.
5.1.

Optimal Power Flow Problem

5.1.1. Description. This problem is concernedwith the management of power transmission over a
network that has observable parts, where measured
data are available, and unobservable parts (Dancre
et al. 2002). For an unobservable part, one usually has
estimations of the data, for example, from the past.
For a proper management of the network, however, acomplete and actual database for the entire networkis needed. Therefore, one relies on computed data in allparts of the network. In the observable areas, the computed data should be as close to the measured dataas possible. In addition, one needs to minimize theweighted least-squares distance between the computed
data and the estimated data in unobservable parts andboundary areas. This gives a nonlinear optimization
problem of the form
minn f x
x∈

s.t.

gx = 0l ≤ hx ≤ u

(9)with the objective function f  n → , the equality
constraints g n → m , and the inequality constraintsh n → p being twice continuously differentiable.
The optimization problem (9) can be solved using
an interior-point based tool such as LOQO or IPOPT(Vanderbei and Shanno 1999, Wächter and Biegler
2006). These solvers require the provision of theJacobian matrices $gx and $hx as well as the
Hessian of the Lagrange function Lx 8 9 = f x +8T gx + 9T hx with respect to x in sparse format.
We report on runtimes for this Hessian computation
using actual problem instances.

5.1.2. Results and Discussion. Table 1 lists theabsolute runtimes in seconds spent in the various
steps for the ﬁve Hessians considered in our experiments. The ﬁrst two columns of Table 2 show thenumber of rows and the average number of nonzeros
per row in the Hessians used in the experiments; thenext two columns show the number of colors used in
the direct and indirect cases; and the last four columns
show timing results for various steps, each normalized relative to the time needed to evaluate the underlying function f .
The results in Tables 1 and 2 clearly show that
employing coloring in Hessian computation enables
one to solve large-size problems that could not otherwise have been solved. For problem sizes where
dense computation is possible, the results show thatsparsity exploitation via coloring yields huge savingsin runtime. Furthermore, it can be seen that indirect computation using acyclic coloring is faster thandirect computation using star coloring, considering
overall runtime. Comparing the steps S1, S2, and S3against each other, as can be seen from Figure 6, thecoloring (S1) and recovery (S3) steps are almost negligible compared with the step in which the Hessianseed matrix product is computed (S2), both in the
direct and indirect methods.
Numerically, we observed that indirect recovery
gave Hessian entries of the same accuracy as direct
recovery. This experimental observation agrees wellwith the analysis in §3.3.
Table 2Matrix Structural Data, Number of Colors, and Normalized
Runtime Relative to Function Evaluation for Test Hessians in
the Optimal Power Flow Problem
ColorsTotal (S1–S3)n

DirectIndirect

S0

DirectIndirect

Dense

72
3,760
4,472
9,932
22,540

3.72
4.11
3.99
3.91
3.94

9
15
16
22
16

6
12
10
9
10

125
421
519
3213
10700

954
1419
1304
1520
9630

725
1285
976
872
7020

10000342964359819680512Note. The asteriks indicate that space could not be allocated.

∗∗∗Gebremedhin et al.: Efﬁcient Computation of Sparse Hessians Using Coloring and AD

219

INFORMS Journal on Computing 21(2), pp. 209–223, © 2009 INFORMS

Indirect method

Direct method
10,000
Total
S2
S3
S1

8,000

Runtime (task)/runtime (F)Runtime (task)/runtime (F)10,000

6,000

4,000

2,000
0

8,000

6,000

4,000

2,000
0

0

5

10

20

0

5

10

n/1,000
Figure 6

Runtimes of the Various Steps Normalized by the Runtime of Function Evaluation for the Power Flow Problem

A ﬁnal point to be noted from Table 2 is that theruntime of the sparsity detection routine is relatively
large in comparison with the routines in the othersteps. In future work, we plan to explore ways inwhich this can be reduced.
5.2.

20

n/1,000

Unconstrained Quadratic Optimization
Problem

5.2.1. Description. The sizes and structures of the
Hessians from the optimal power ﬂow problem that
we could include in our experiments were quite limited. To be able to study the performance of thevarious steps in a systematic fashion, we considereda synthetic problem in which we have a direct control
over the size and structure of the Hessians. In particular, we used the unconstrained quadratic optimization
problem minx∈n f x with f x = xT Cx +aT x, C ∈ n×n
and aT = 10     10 ∈ n , where the Hessian is simply the matrix C. We considered two kinds of sparsity
structures for the matrix C: banded (bd) and random(rd). Furthermore, the test matrices were designed insuch a way that:(i) the number of nonzeros per row in a bandedmatrix (denoted by :) is nearly the same as the number of nonzeros per row in a random matrix (denotedby :),
 and(ii) the value for :, or :, remains constant as the
problem dimension n is varied.
In our experiments, we used the values : ∈10 20, : ∈ 1098 2099, and n/1000 ∈  ≡5 10 20 40 60 80 100.
5.2.2. Results and Discussion.
Number of Colors. Table 3 provides a summary of
the numbers of colors used by the star-coloring algorithm (direct method) and the acyclic-coloring algorithm (indirect method) for all the sparsity structuresand input sizes considered in the experiments. Two
observations can be made from this table.
First, for the banded structure, the acyclic- andstar-coloring algorithms invariably used :/2 + 1 and

2:/2 + 1 colors, respectively, regardless of the valueof n. In view of Theorems 4.3 and 4.4, and noting that
:/2 is the bandwidth of the corresponding graphs,
we see that both algorithms ﬁnd optimal solutions forband graphs. Both algorithms are greedy, and vertices
were colored in the natural ordering of the graphs.
Hence, the observed phenomenon agrees with thetheory of distance-1 coloring discussed in the lastparagraph of §4.
Second, in both the star- and the acyclic-coloringcases, the numbers of colors required by the random
structures were observed to be nearly twice the corresponding numbers in the banded structures. Moreover, the numbers of colors varied only slightly as the
problem dimension n was varied.
Runtime. Table 4 lists the absolute runtime in seconds spent in the various steps while using a direct
and an indirect method. The information in Table 4
is analogous to that presented in Table 1 for the optimal power ﬂow problem. The general conclusion to
be drawn from Table 4 in terms of the enabling powerof the coloring techniques in the overall computation
is similar to that drawn from the optimal power ﬂow
problem. Our objective here is to show how the execution time for each step grows as a function of the
input size.
Figure 7 shows a collection of normalized runtime
versus problem dimension (n) plots. In particular, theTable 3

Number of Colors Used by the Star- andAcyclic-Coloring Algorithms for All
Problem Dimensions n in the Set n/1000 ∈ ≡ 5 10 20 40 60 80 100
, 10, 10.98

bd
rd20, 20.99

Star

Acyclic

Star

Acyclic

11
21–24

6
9–11

21
50–56

11
18–19

Gebremedhin et al.: Efﬁcient Computation of Sparse Hessians Using Coloring and AD

220

INFORMS Journal on Computing 21(2), pp. 209–223, © 2009 INFORMS

Table 4Absolute Runtimes in Seconds for the Evaluation of the Function f and the Steps S1, S2, and S3
in the Quadratic Optimization Problem
DirectIndirectn/1000

Evalf 

S1

S2

S3

Tot

S1

S2

S3

Tot

Dense

bd:
5
10
20
40
60
80
100

0.00060.00100.00190.00350.00490.00620.0077

005
011
023
055
081
103
131

020
040
079
161
231
306
383

0.006
0.012
0.023
0.046
0.066
0.088
0.110

025
052
105
221
319
418
525

008
025
081
286
612
1103
1633

011
022
043
086
125
165
208

005
010
021
046
073
099
123

024
057
145
418
810
1367
1964

861
3420

rd:
5
10
20
40
60
80
100

0.00050.00090.00200.00430.00820.01250.0170

007
019
046
108
177
250
326

041
090
201
451
754
1078
1316

0.007
0.014
0.029
0.064
0.103
0.140
0.174

049
110
250
566
941
1341
1660

010
030
095
321
675
1152
1752

017
038
083
184
313
454
555

011
030
087
267
525
866
1301

038
098
265
773
1513
2471
3607

1479
5894

bd:
5
10
20
40
60
80
100

0.00080.00150.00290.00560.00840.01090.0135

016
036
084
168
256
354
444

063
127
257
514
746
1047
1295

0.015
0.030
0.061
0.118
0.176
0.234
0.300

081
166
347
693
1019
1424
1770

014
037
107
340
693
1185
1774

033
067
134
261
391
547
681

012
023
051
110
166
231
307

059
128
291
712
1251
1963
2761

840
3628

rd:
5
10
20
40
60
80
100

0.00080.00150.00330.00870.01500.0246
0.0317

029
070
165
391
604
891
1193

166
349
814
1999
2991
4512
5471

0.017
0.035
0.071
0.156
0.234
0.338
0.413

197
423
986
2406
3619
5436
6707

020
052
144
435
856
1437
2134

055
115
267
643
989
1500
1900

024
065
187
559
1081
1803
2691

099
231
598
1636
2926
4740
6725

1510
5949

∗∗∗
∗∗∗
∗∗∗
∗∗∗
∗∗∗

∗∗∗
∗∗∗
∗∗∗
∗∗∗
∗∗∗

∗∗∗
∗∗∗
∗∗∗
∗∗∗
∗∗∗

∗∗∗
∗∗∗
∗∗∗
∗∗∗
∗∗∗

Notes. The upper half of the table shows results for  = 10 (banded) and  = 1098 (random) and the lower halffor  = 20, and  = 2099. All runtimes are averages of ﬁve runs. For the random structures, the runtimes are in
addition averaged over ﬁve randomly generated matrices.

vertical axis in each subﬁgure shows the execution
time of a speciﬁc step divided by the time needed toevaluate the function f being differentiated; note thatthe scales on the axes differ from subﬁgure to subﬁgure. Below, we discuss the runtime behavior of thevarious steps turn by turn; but ﬁrst, we look at how
the normalizing quantity, the time needed for evaluating f , itself grows as a function of n.
Time for Evaluating f . Because the number ofnonzeros per row (column) in the structures we considered is constant, the time needed to evaluate thefunction f theoretically is expected to be linear in
the number of rows (columns) n. Figure 8 shows thatthe practically observed execution times are roughly
linear in n across the structures we considered. For
the banded structures, the growth is actually slightly
sublinear. The growth is somewhat superlinear forthe random structures, especially for the cases where: = 2099. This is due mainly to the irregular memory

accesses involved and the associated nonuniform
costs in hierarchical memory.
Step S1: Coloring and Generation of Seed Matrix.
Recall from §2 that the complexity of the star-coloringalgorithm for a graph on n vertices is Ond¯2 andthat of the acyclic-coloring algorithm is Ond¯2 · ,where is the inverse of Ackermann’s function. For
the banded sparsity structures, the quantity d¯2 inthe associated adjacency graphs is nearly :2 , independent of the parameter n. In light of these facts,
the trends observed in the various cases in Figure 7
are in agreement with theoretical analyses. For the
banded structures (the top two rows), it can be seenthat the runtime of the star-coloring algorithm grows
linearly with n (left column), whereas the runtimeof the acyclic-coloring algorithm is slightly superlinear (right column). The general trend in the random
structures is very similar, but slightly more erratic,
again due to irregular memory accesses.

Gebremedhin et al.: Efﬁcient Computation of Sparse Hessians Using Coloring and AD

221

INFORMS Journal on Computing 21(2), pp. 209–223, © 2009 INFORMSbd, ρ = 10, direct method

bd, ρ = 10, indirect method
2,500

Total
S2
S1
S3

2,000

Runtime (task)/runtime (F)Runtime (task)/runtime (F)

2,500

1,500
1,000
500
0

5 10

20

40

60

80

500

5 10

20

40

bd, ρ = 20, direct method

bd, ρ = 20, indirect method

1,000

500

5 10

20

40

60

80

80

100

80

100

80

100

1,000

500

5 10

20

40

60

n/1,000

n/1,000rd, ρ = 10.98, direct methodrd, ρ = 10.98, indirect method

Runtime (task)/runtime (F )

2,000

1,500

1,000

500

5 10

20

40

60

80

1,500

1,000

500

0

100

5 10

20

40

n/1,000

60

n/1,000rd, ρ = 20.99, direct methodrd, ρ = 20.99, indirect method

4,000

4,000

Runtime (task)/runtime (F)Runtime (task)/runtime (F)

100

1,500

0

100

0

3,000

2,000

1,000

5 10

20

40

60

n/1,000
Figure 7

80

2,000

1,500

0

60

n/1,000

Runtime (task)/runtime (F )Runtime (task)/runtime (F )1,000

n/1,000

2,000

Runtime (task)/runtime (F )

1,500

0

100

2,000

0

2,000

80

100

3,000

2,000

1,000

0

5 10

20

40

60

n/1,000

Execution Time of the Various Steps Normalized by the Time Needed for Function Evaluation vs. Problem Size

Gebremedhin et al.: Efﬁcient Computation of Sparse Hessians Using Coloring and AD

222

INFORMS Journal on Computing 21(2), pp. 209–223, © 2009 INFORMS× 107
bd, ρ = 10bd, ρ = 20rd, ρ = 10.98
rd, ρ = 20.98

3.5

Runtime (f )/n

3.0
2.5
2.0
1.5
1.0
0.5

5 10

20

40

60

80

100

n/1,000
Figure 8

Runtime of Function Evaluation in Seconds Normalized by
Input Size n vs. n

Step S2: Computation of the Compressed
Hessian. Figure 7 shows that the time for the stepin which the compressed Hessian HS is computed
is linear in the problem dimension n. The analyticaljustiﬁcation for this behavior stems from two sources.
First, as mentioned earlier, the time needed for evaluating the function f is linear in n. Second, the numberof columns p in a seed matrix (the number of colorsused) remained constant or nearly constant as the
problem dimension n in our experiments was varied,both in a direct and an indirect method. Theoretically,
the complexity of computing the Hessian-seed vector
product using AD is known to be a small constant(in the order of 10) times the time need to evaluate
the function being differentiated (Griewank 2000).
Hence, the fact that the observed runtime grewlinearly with n for both structures is consistent withtheoretical analyses.
Step S3: Recovery of the Original Hessian Entries.
As discussed in §3, the complexity of directRecover1 is Om: , where m is the number of nonzeros in the Hessian. The constant hidden in thisexpression is rather small, because the computation
involved is fairly easy. Similarly, the complexity of
indirectRecover, which relies on the use of twocolored trees, was shown to be Om . Due to the overhead associated with the management of nontrivial
data structures, the hidden constant here is expected
to be considerably larger, to the extent that the execution time of the routine in practice becomes morethan the corresponding time for directRecover1. The
observed runtimes in Figure 7 clearly reﬂect thesefacts. For a similar reason, even though directRecover2 theoretically is faster than directRecover1,
we used the latter in our experiments because it islikely to be faster in practice.
Overall Runtime. Considering all the steps together, is a direct method faster or slower than an

indirect method? The results in Figure 7 show thatthe answer depends on the size and structure of
the Hessian being computed. For the random structure with nearly 20 nonzeros per row, an indirect
method is consistently observed to be faster than adirect method. A similar statement can be made forthe banded structures of relatively small size (n upto 20000). For larger-size banded problems and formany of the random matrices with nearly 10 nonzeros
per row, a direct method was observed to be faster.
These observations are in contrast to those in the optimal power ﬂow problem where an indirect method
was always found to be faster. Comparing the relative
contribution of the various steps with the total runtime, we observe that the Hessian-seed matrix product step (S2) is by a large margin the most expensive
in a direct method, whereas the coloring step (S1) isslightly the dominant step in an indirect method.
Numerical Accuracy. As in the optimal power
ﬂow problem, here again, the numerical values ofthe Hessian entries obtained using indirectRecoverwere observed to be of the same accuracy as the values obtained using directRecover1—a typical pairof values obtained using the two methods matched in
all of the computed digits in double precision.

6.

Conclusion

We studied compression-based calculation of sparse
Hessians using automatic differentiation. We considered the case where a matrix is compressed suchthat the recovery is direct (star coloring) and thecase where the recovery requires additional arithmetic
work (acyclic coloring). Our experimental results
showed that sparsity exploitation via star and acyclic
coloring enables one affordably to compute Hessiansof dimensions that could not have been computed
otherwise. For sizes where a computation that does
not exploit sparsity is at least possible, the results
showed that the techniques render dramatic savingsin runtime. We believe savings of similar magnitude
would be attained should an AD tool other than
ADOL-C be used, because the execution time for
the Hessian-seed matrix product is likely to dominate the overall runtime for any reasonable function.
The experimental results also showed that, for realworld optimization problems, an acyclic-coloringbased method is faster than a star-coloring-based
method, considering the overall process. Furthermore,
we showed, both analytically and experimentally, thatindirect recovery using two-colored trees is numerically as stable as direct recovery.
Acknowledgments

The authors thank Dr. Fabrice Zaoui of EDF R&D MOSSE,
Clamart, France for helping with the experiments on the
power ﬂow problem. They also thank the anonymousGebremedhin et al.: Efﬁcient Computation of Sparse Hessians Using Coloring and AD
INFORMS Journal on Computing 21(2), pp. 209–223, © 2009 INFORMS

referees for their valuable comments, which helped them
improve the quality of the paper. This work was supportedby the Ofﬁce of Science of the U.S. Department of Energy
under the Scientiﬁc Discovery Through Advanced Computing (SciDAC) program through Grant DE-FC-0206-ER-25774
awarded to the CSCAPES Institute, by the U.S. National Science Foundation Grant ACI 0203722, and by the German
Research Foundation Grant Wa 1607/2-1. A. H. Gebremedhin’s current afﬁliation is Department of Computer Science and Computing Research Institute, Purdue University,
West Lafayette, IN 47907.

References
Brandstädt, A., V. B. Le, J. P. Spinrad. 1999. Graph Classes: A Survey: Monographs on Discrete Mathematics and Applications. SIAM,
Philadelphia.
Büskens, C., H. Maurer. 2001. Sensitivity analysis and real-time
optimization of parametric nonlinear programming problems.
M. Gröschel, S. Krumke, J. Rambau, eds. Online Optimization of
Large Scale Systems. Springer, Berlin, 3–16.
Coleman, T. F., J. Cai. 1986. The cyclic coloring problem and estimation of sparse Hessian matrices. SIAM J. Algebraic Discrete
Methods 7(2) 221–235.
Coleman, T. F., J. J. Moré. 1984. Estimation of sparse Hessian
matrices and Graph coloring problems. Math. Programming 28
243–270.
Dancre, M., P. Tournebise, P. Panciatici, F. Zaoui. 2002. Optimal
power ﬂow applied to state estimation enhancement. Proc. 14th
Power Systems Comput. Conf., Sevilla, Spain, Paper 3, Session
37.
Dixon, L. 1991. Use of automatic differentiation for calculating Hessians and Newton steps. A. Griewank, G. Corliss, eds. Automatic Differentiation of Algorithms: Theory, Implementation, andApplication, SIAM, Philadelphia, 114–125.
Gay, D. 1996. More AD of nonlinear AMPL models: Computing Hessian information and exploiting partial separability. M. Berz,

223

C. Bischof, G. Corliss, A. Griewank, eds. Computational Differentiation: Techniques, Applications, and Tools. SIAM, Philadelphia,
173–184.
Gebremedhin, A. H., F. Manne, A. Pothen. 2005. What color is your
Jacobian? Graph coloring for computing derivatives. SIAM Rev.
47(4) 629–705.
Gebremedhin, A. H., A. Tarafdar, D. Nguyen, A. Pothen. 2009.
ColPack: A graph coloring package for supporting sparse derivative matrix computation. ACM Trans. Math. Software. In
preparation.
Gebremedhin, A. H., A. Tarafdar, F. Manne, A. Pothen. 2007. New
acyclic and star coloring algorithms with application to computing Hessians. SIAM J. Sci. Comput. 29 1042–1072.
Griewank, A. 2000. Evaluating Derivatives: Principles and Techniquesof Algorithmic Differentiation. Frontiers in Applied Mathematics,
Number 19. SIAM, Philadelphia.
Griewank, A., D. Juedes, J. Utke. 1996. ADOL-C: A package forthe automatic differentiation of algorithms written in C/C++.
ACM Trans. Math. Software 22 131–167.
Lin, Y., S. S. Skiena. 1995. Algorithms for square roots of graphs.
SIAM J. Discrete Math. 8 99–118.
Lovász, L. 1972. A characterization of perfect graphs. J. Combin.
Theory 13 95–98.
McCormick, S. T. 1983. Optimal approximation of sparse Hessians
and its equivalence to a graph coloring problem. Math. Programming 26 153–171.
Papadimitriou, C. 1976. The NP-completeness of the bandwidth
minimization prolem. Computing 16 263–270.
Powell, M. J. D., P. L. Toint. 1979. On the estimation of sparse
Hessian matrices. SIAM J. Numer. Anal. 16(6) 1060–1074.
Vanderbei, R., D. Shanno. 1999. An interior-point algorithm fornonconvex nonlinear programming. Comput. Optim. Appl. 13
231–252.
Wächter, A., L. Biegler. 2006. On the implementation of an interiorpoint ﬁlter line-search algorithm for large-scale nonlinear programming. Math. Programming 106(1) 25–57.
Walther, A. 2008. Computing sparse Hessians with automatic differentiation. ACM Trans. Math. Software 34(1) Paper 3.