The Efficient Computation of Sparse Jacobian Matrices Using Automatic Differentiation.
This paper is concerned with the efficient computation of sparse Jacobian matrices of nonlinear vector maps using automatic differentiation (AD). Specifically, we propose the use of a graph coloring technique, bicoloring, to exploit the sparsity of the Jacobian matrix J and thereby allow for the efficient determination of J using AD software. We analyze both a direct scheme and a substitution process. We discuss the results of numerical experiments indicating significant practical potential of this approach.
Introduction
The efficient numerical solution of nonlinear systems of algebraic equations, F
usually requires the repeated calculation or estimation of the matrix of first derivatives, the Jacobian
In large-scale problems matrix J is often sparse and it is important to exploit
this fact in order to efficiently determine, or estimate, matrix J at a given argument x. This paper
is concerned with the efficient calculation of sparse Jacobian matrices by the judicious application of
automatic differentiation techniques. Specifically, we show how to define "thin" matrices V and W
such that the nonzero elements of J can easily be extracted from the calculated pair (W T J; JV ).
This research was partially supported by the Applied Mathematical Sciences Research Program (KC-04-02) of
the Office of Energy Research of the U.S. Department of Energy under grant DE-FG02-90ER25013, and in part by
the Advanced Computing Research Institute, a unit of the Cornell Theory Center which receives major funding from
the National Science Foundation and IBM Corporation, with additional support from New York State and members
of its Corporate Research Institute. This technical report also appears as Cornell Computer Science Technical Report
TR 95-1557.
y Computer Science Department and Center for Applied Mathematics, Cornell University, Ithaca NY 14850.
z Computer Science Department, Cornell University, Ithaca NY 14850.
Given an arbitrary n-by-t V matrix V , product JV can be directly calculated using automatic
differentiation in the "forward mode"; given an arbitrary m-by-t W matrix W , the product W T J
can be calculated using automatic differentiation in the "reverse mode", e.g., [11, 13].
The forward mode of automatic differentiation allows for the computation of product JV in
time proportional to t V \Delta !(F ) where !(F ) is the time required to evaluate F . This fact leads to
the following practical question. Given the structure of a sparse Jacobian matrix J , how can a
matrix V be chosen so that the nonzeros of J can easily be determined from the product JV ? A
good solution is offerred by the sparse finite-differencing literature [4, 5, 6, 7, 8, 10] and adapted to
the automatic differentiation setting [1]. Partition the columns of J into a set of groups GC , where
the number of groups in GC is denoted by jGC j, such that the columns in each group G 2 GC are
structurally orthogonal 1 . Each group G 2 GC defines a column v of only if column
i is in group G; otherwise, v It is clear that the nonzeros of J can be immediately "identified"
from the computed product JV . Graph coloring techniques, applied to the column intersection
graph of J , can be used to try and produce a partition GC with low cardinality jGC j. This, in turn,
induces a thin matrix V , i.e., construct However, it is not always
possible to ensure that jGC j is small: consider a sparse matrix J with a single dense row.
Alternatively, the reverse mode of automatic differentiation allows for the computation of the
product W T J in time proportional to t W \Delta !(F ) where !(F ) is the time required to evaluate F , and
. The "transpose" of the argument above can lead to an efficient way to determine J .
That is, apply graph coloring techniques to the row intersection graph of J to induce a thin matrix
via the reverse mode of automatic differentiation (takes time proportional to
trivially extract the nonzeros of J from the computed matrix W T J . Of course it is easy
to construct examples where defining a thin matrix W is not possible - consider the case where J
has a dense column.
Clearly there are problems where a row-oriented approach is preferable, there are problems
where a column-oriented approach is better. Unfortunately, it is easy to devise problems where
neither approach is satisfactory: let J have both a dense row and a dense column. This is exactly
when it may pay to use both modes of automatic differentiation simultaneously: compute a pair
choices of W and V , and extract the nonzero elements of J from this
computed pair of thin matrices.
Our concern in this paper is with efficiency with respect to number of floating point operations
or f lops. We do not concern ourselves with space requirements in this study. However, it should
be noted that the reverse mode of automatic differentiation often requires significantly more space
than the forward mode: if space is tight then our suggested approach, which involves application
of both forward and reverse modes, may not be possible. There is current research activity on
reducing the space requirements of the reverse mode of automatic differentiation, e.g., [12].
We note that an independent proposal regarding sparse Jacobian calculation is made by Hossain
and Steihaug [15]: a graph-theoretic interpretation of the direct determination problem is given and
an algorithm based on this interpretation is provided. In this paper we proffer a new direct method
and we also propose a substitution method, both based directly on the Jacobian structure. We
compare our direct and substitution methods, numerically, and we discuss the round-off properties
of the substitution method. In addition, we interface our graph coloring software to the automatic
differentiator ADOL-C, [14], and report on a few preliminary computational results.
The remainder of the paper is organized as follows. In x2, we review the relevant aspects of
Two nonzero n-vectors v; w are structurally orthogonal if v i   w
automatic differentiation, both forward and reverse modes. In x3, we formalize the combinatorial
problems to be solved both from a matrix point of view and in terms of graph theory. We propose
both a direct determination problem and a substitution problem. In x4, we propose "bi-coloring"
approaches to both the direct determination and "determination by substitution" problems. The
bi-coloring technique produces matrices V and W where JV is subsequently determined via the
forward mode of automatic differentiation, W T J is detemined via the reverse mode. Typically, the
column dimensions of V and W will be small: the cost of the application of automatic differentiation
is proportional to the sum of the column dimensions of V and W (times the work to evaluate F ).
In x5 we present and discuss various numerical experiments. The experiments indicate that our
bi-coloring approach can significantly reduce the cost of determining J (over one-sided Jacobian
determination).
The substitution method we propose consistently outperforms the direct method. However, the
substitution calculation increases the chance of round-off contamination. This effect is discussed in
x6. We end the paper in x7 with some concluding remarks and observations on possible directions for
future research. Specifically, we note that while sparsity is a symptom of underlying structure in a
nonlinear problem, it is not a necessary symptom. Moreover, it is often possible to exploit structure
in the absence of sparsity and apply AD tools "surgically" to efficiently obtain the Jacobian matrix
J .
/sectionBasics of automatic differentiation Automatic differentiation is a chain rule based technique
for evaluating the derivatives analytically (and hence without any truncation errors) with
respect to input variables of functions defined by a high-level language computer program. In this
section we briefly review the basics of automatic differentiation, borrowing heavily from [11, 13].
A program computing the function can be viewed as a sequence
of scalar assignments v where the vector v can be thought of as set of ordered
variables such that v i computed before v j using the set of variables fv k jk ! ig. Here / j
represent elementary functions, which can be arithmetic operations and/or univariate transcedental
functions.
Ordering the variables as above, we can partition variables v j into three vectors.
In general, the number of intermediate variables is much larger than the dimensions of the problem,
Assume that all these elementary functions / i are well defined and have continuous elementary
Assuming with out loss of generality that the dependent variables y do
not themselves occur as arguments of elementary functions, we can combine the partials c ij into
the (p +m) \Theta (n + p) matrix
c n+i;j
1-i-p+m
Unless elementary functions with more than two arguments are included in the library, each
row of C contains either one or two non zero entries. We define a number
Also, since work involved in an elementary function is proportional to number of arguments, it
follows that
proportionality. Because of the ordering relation the square matrix C is upper
trapezoidal with nontrivial superdiagonals. Thus C can be partitioned as
A L
where
and L is strictly lower triangular. Application of the chain rule yields:
!/
\Deltax
\Deltay
If we eliminate the intermediate vector \Deltay from above (1.4), we get an expression(the Schur
complement) for the Jacobian:
is a unit lower triangular matrix, the calculation of the matrix products ~
leads to two natural ways to compute J :
A or
The alternative expressions for J given in (1.6) define the two basic modes of automatic differenti-
ation, forward and reverse.
The forward mode corresponds to computing the rows of ~
A, one by one, as the corresponding
rows of [A; are obtained from successive evaluation of elementary functions. Since this
amounts to solutions of n linear systems with lower-triangular matrix [I \Gamma L], followed by multiplication
of dense columns of ~
A by M , the total computational effort is roughly n \Delta q or n \Delta !(F ).
The reverse mode corresponds to computing ~
M as solution to linear system
M T . This back-substitution process can begin only after all elementary functions and their partial
derivatives have been evaluated. Since this amounts to the solution of m linear systems with lower
triangular matrix [I \Gamma L], followed by multiplication of dense rows of ~
M by A, total computational
effort is roughly m \Delta q or m \Delta !(F ).
We are interested in computing products of the form JV and W T J . Product JV can be
computed:
which can clearly be done in time proportional to t V n\Thetat V . Analagously, product
can be computed:
which can be done in time proportional to t W \Delta !(F ) assuming
An important subcase worthy of special attention is when F is a scalar function, i.e.,
In this case the Jacobian matrix corresponds to the transpose of the gradient of F and is a single
row vector. Note that the complexity arguments applied to this case imply that the reverse mode
of AD yields the gradient in time proportional to !(F ) whereas the forward mode costs n \Delta !(F ).
The efficiency of the forward mode evaluation of the gradient can be dramatically increased - i.e.,
the dependence on n is removed - if F has structure that can be exploited [2].
Partition problems and graph theory
Our basic task is to efficiently determine thin matrices V; W so that the nonzero elements of J can
be readily extracted from the information (W T J; JV ). The pair of matrices (W T J; JV ) is obtained
from the application of both modes of automatic differentiation: matrix W T J is computed by the
reverse mode, the forward mode determines JV . The purpose of this section is to more rigorously
formulate the question of determining suitable matrices V; W , first in the language of "partitions",
and then using graph theoretic concepts.
We begin with an example illustrating the usefulness of simultaneously applying both modes of
AD, forward and reverse. Consider the following n-by-n Jacobian, symmetric in structure but not
in value:
It is clear that a partition of columns consistent with the direct determination of J requires n
groups. This is because a "consistent column partition" requires that each group contain columns
that are structurally orthogonal and the presence of a dense row implies each group consists of
exactly one column. Therefore, if matrix V corresponds to a "consistent column partition" then V
has n columns and the work to evaluate JV by the forward mode of AD is proportional to n \Delta !(F ).
By a similar argument, and the fact that a column of J is dense, a "consistent row partition"
requires n groups. Therefore, if matrix W corresponds to a "consistent row partition" then W has
rows and the work to evaluate W T J by the reverse mode of AD is proportional to n \Delta !(F ).
Definition 2.1 A bi-partition of a matrix A is a pair (GR ; GC ) where GR is a row partition of
a subset of the rows of A , GC is a column partition of a subset of the columns of A.
In this example the use of a bi-partition dramatically decreases the amount of work required to
determine J . Specifically, the total amount of work required is proportional to 3 \Delta !(F ). To see this
the usual convention of representing
the i th column of the identity matrix with e i . Clearly elements are directly determined from
the product JV ; elements 4 are directly determined from the product W T J .
The basic idea is to partition the rows into a set of groups GR and the columns into a set of
groups GC , with small as possible, such that every nonzero element of J can be
directly determined from either a group in GR or a group in GC .
Definition 2.2 A bi-partition (GR ; GC ) of a matrix A is consistent with direct determination
if for every nonzero a ij of A, either column j is in a group of GC which has no other column having
a nonzero in row i, or row i is in a group of GR which has no other rows having a nonzero in column
j.
Clearly, given a bi-partition (GR ; GC ) consistent with direct determination, we can trivially
construct matrices n\ThetajG C j such that A can be directly determined from
If we relax the restriction that each nonzero element of J be determined directly then it is
possible that the work required to evaluate the nonzeroes of J can be further reduced. For example
we could allow for a "substitution" process when recovering the nonzeroes of J from the pair

Figures

2.1, 2.2 illustrate that a substitution method can win over direct determination:

Figure

2.1 corresponds to direct determination, Figure 2.2 corresponds to determination using

Figure

2.1: Optimal partition for direct method
In both cases elements labelled   are computed from the column grouping, i.e., calculated using
the product JV ; elements labelled 4 are calculated form the row groupings, i.e., calculated using
the product W T J . The matrix in Figure 2.1 indicates that we can choose GC with jGC
GR with jGR determine all elements directly. That is, choose choose
Therefore in this case the work to compute J satisfies !(J) - 3 \Delta !(F ). Note
that some elements can be determined twice, e.g., J 11 .
However, the matrix in Figure 2.2 shows how to obtain the nonzeroes of J , using substitution,
in work proportional to 2 \Delta !(F ). Let be the (forward) computed

Figure

2.2: Optimal partition for substitution method
vector
2 be the (reverse) computed row vector p T
Most of the nonzero elements are determined directly (no conflict). The remaining elements can
be resolved,
It is easy to extend this example so that the difference between the number of groups needed,
between substitution and direct determination, increases with the dimension of the matrix. For
example, a block generalization is illustrated in Figure 2.3: if we assume l ? 2w it is straightforward
to verify that in the optimal partition the number of groups needed for direct determination will
be 3w and determination by substitution requires 2w groups.
Definition 2.3 A bi-partition (GR ; GC ) of a matrix A is consistent with determination by
substitution, if there exists an ordering - on elements a ij , such that for every nonzero a ij of A,
either column j is in a group where all nonzeros in row i, from other columns in the group, are
ordered lower than a ij , or row i is in a group where all the nonzeros in column j, from other rows
in the group, are ordered lower than a ij .
In the usual way we can construct a matrix V from the column grouping GC and a matrix W
from the row grouping GR : for example, to construct the columns of V associate with each group
l
l

Figure

2.3: Block example
in GC a boolean vector, with unit entries indicating membership of the corresponding columns. We
can now state our main problem(s) more precisely:
The bi-partition problem (direct) : Given a matrix A, obtain a bi-partition (GR ; GC )
consistent with direct determination, such that total number of groups,
minimized.
The bi-partition problem (substitution) : Given a matrix A, obtain a bi-partition (GR ; GC )
consistent with determination by substitution, such that total number of groups,
jGC j, is minimized.
The bi-partition problems can also be expressed in terms of graphs and graph coloring. This
graph view is important in that it more readily exposes the relationship of the bi-partition problems
with the combinatorial approaches used in the sparse finite-differencing literature, e.g., [4, 5, 6, 7, 8].
However, we note that the remainder of this paper, with the exception of the error analysis in x6,
does not rely directly on this graph interpretation.
To begin, we need the usual notion of a coloring of the vertices of a graph, the definition of a
bipartite graph, and the concept of path coloring [4, 7] specialized to the bipartite graph case.
A p-coloring of a graph is the set of vertices or nodes and E is the set of
edges, is a function
such that OE(u) 6= OE(v), if (u; v) 2 E. The chromatic number -(G) is the smallest p for which G has
a p-coloring. A p-coloring OE of G induces a partition of vertices
such that
Given a matrix A 2 ! m\Thetan , define a bipartite graph E) where
corresponds to the jth column of A, and r i corresponds to
the ith row of A. There is an edge, (r a nonzero in A.
In [4, 7] a path p-coloring of a graph is defined to be a vertex coloring using p colors with
the additional property that every path of at least 3 edges uses at least 3 colors. Here we need a
slight modification of that concept appropriate for the direct determination problem. We note that
"color 0" is distinguished in that it corresponds to the lack of a true color assignment: i.e.,
indicates that vertex i is not assigned a color.
Definition 2.4 Let G E) be a bipartite graph. A mapping OE
is a bipartite path p-coloring of G b if
1. Adjacent vertices have different assignments, i.e., if (i;
2. The set of positive colors used by vertices in V 1 is disjoint from the set of positive colors
used by vertices in 0g.
3. If vertices i and j are adjacent to vertex k with
4. Every path of 3 edges uses at least at least 3 colors.
The smallest number for which graph G b is bipartite path p-colorable is denoted by -

Figure

2.4 shows a valid bipartite path p-coloring. Numbers adjacent to the vertices denote colors.
We note that Hossain and Steihaug [15] define a similar concept. However, their definition of path
p-coloring does not allow for the "uncolor assignment", i.e., Consequently, a technique to
remove empty groups is needed [15].
We are now in position to state the graph analogy to the concept of a bi-partition consistent
with direct determination.
Theorem 2.1 Let A be a m \Theta n matrix with corresponding bipartite graph E).
The mapping induces a bi-partition (GR ; GC ), with
consistent with direct determination if and only if OE is a bipartite path p-coloring of G b (A).
Proof. (() Assume that OE is a bipartite path p-coloring of G b (A), inducing a bi-partition (GR ; GC )
of rows and columns of A. If this bi-partition is not consistent with direct determination, then there
is a nonzero element a ij in the matrix for which the definition "either column j is in a group of
GC which has no other column having a nonzero in row i, or row i is in a group of GR which has
no other rows having a nonzero in column j" doesn't hold. This can happen only if one of the
following cases hold:
there exists a column q with a iq 6= 0, such that OE(c j
this contradicts condition 3.
there exists a row p with a pj 6= 0, such that OE(r i
contradicts condition 3.
COLUMNS
ROWS000000

Figure

2.4: A valid bipartite path coloring
There exists a column q and a row p, such that columns j and q are
in the same group with a iq 6= 0 and rows i and p are in the same group with a pj 6= 0 This
implies OE is a 2-coloring of the path (r contradiction of condition 4.
Conversely, assume that OE induces a bi-partition consistent with direct determination of
A. It is clear that conditions 1 \Gamma 3 must be satisfied. It remains for us to establish condition 4: i.e.,
every path of 3 edges uses at least 3 colors. Suppose there is a bi-colored
by condition 3 the two colors on this path are positive.
It is easy to see that element a jk cannot be determined directly: there is a conflict in row group
there is a conflict in column group OE(c k these are the only two
chances to determine a jk .  .
To capture the substitution notion the cyclic p-coloring definition [4] is modified slightly and
applied to a bipartite graph.
Definition 2.5 Let A be a m \Theta n matrix with corresponding bipartite graph E).
A mapping OE is a bipartite cyclic p-coloring of G b if
1. Adjacent vertices have different assignments, i.e., if (i;
2. The set of positive colors used by vertices in V 1 is disjoint from the set of positive colors
used by vertices in 0g.
3. If vertices i and j are adjacent to vertex k with
4. Every cycle uses at least at least 3 colors.
The smallest number for which graph G b is bipartite cyclic p-colorable is denoted by - c (G b ).

Figure

2.5 shows a valid bipartite cyclic p-coloring note that only 2 colors are necessary whereas
the bipartite path p-coloring in Figure 2.4 requires 3 colors. The notion of a bi-partition consistent12COLUMNS

Figure

2.5: A valid bipartite cyclic coloring
with determination via substitution can now be cleanly started in graph-theoretic terms.
Theorem 2.2 Let A be a m \Theta n matrix with corresponding bipartite graph E).
The mapping OE : induces a bi-partition (GR ; GC ), with
consistent with determination by substitution if and only if OE is a bipartite cyclic p-coloring of
G b (A).
Proof. ()) Assume OE induces a bi-partition consistent with determination by substitution but
OE is not a bipartite cyclic p-coloring of G b (A). Clearly condition must hold; it is easy to
see that if condition 3 doesn't hold then not all nonzero elements can be determined. The only
non-trivial violation is condition 4: there is a cycle which has only two colors, i.e all the vertices
in the cycle have the same color c 1 , and all the vertices in the cycle have the same color
c 2 . Note that neither c 1 and c 2 can be equal to 0, since a node colored 0 in a cycle would imply
that its adjacent vertices are both colored differently, implying that there are at least 3 colors.
Consider the submatrix A s of A, corresponding to this cycle. Submatrix A s has at least two non
zeros in each row and in each column, since each vertex has degree 2 in the cycle. But since we are
considering substitution methods only, at least one element of A s needs to be computed directly.
Clearly there is no way to get any element of this submatrix directly, a contradiction.
Conversely, assume that OE is a bipartite cyclic p-coloring of G b (A) but that bi-partition
induced by OE is not consistent with determination by substitution. But, edges (nonzeros)
with one end assigned color "0" can be determined directly: by the definition of bi-coloring there
will be no conflict. Moreover, every pair of positive colors induces a forest (i.e., a collection of
trees); therefore, the edges (nonzeros) in the induced forest can be resolved via substitution [4].
The two bi-partition problems can now be simply stated in terms of optimal bipartite path and
cyclic p-colorings:
The bipartite path p-coloring problem : Determine a bipartite path p-coloring of G b (A)
with the smallest possible value of p, i.e.,
The bipartite cyclic p-coloring problem : Determine a bipartite cyclic p-coloring of G b (A)
with the smallest possible value of p, i.e.,
The graph theoretic view is useful for both analyzing the complexity of the combinatorial
problem and suggesting possible algorithms, exact or heuristic. In fact, using the the p-coloring
notions discussed above, and an approach similar to that taken in [4], it is easy to show that
corresponding decision problems are NP-complete.
Bipartite cyclic p-coloring decision problem (CCDP): Given an integer p - 3 and an arbitrary
bipartite graph G, is it possible to assign a cyclic p-coloring to nodes of G?
Bipartite path p-coloring decision problem (PCDP) : Given an integer p - 3 and an arbitrary
bipartite graph G, is it possible to assign a bipartite path p-coloring to nodes of G?
The proofs are a straightforward adaptation of those in [4] and we omit them here. The upshot
of these (negative) complexity results is that in practise we must turn our attention to (fast)
heuristics to approximately solve the cyclic and path coloring problems. In the next section we
present simple, effective, and "easy-to-visualize" heuristics for these two combinatorial problems.
Finally, it is easy to establish a partial ordering of chromatic numbers:
where G(M) refers to the column intersection graph of matrix M , -(G(M)) is the (usual) chromatic
number of graph G(M ).
The first inequality in (2.2) holds because if OE is a bipartite path p-coloring then OE is a bipartite
cyclic p-coloring; the second inequality holds because a trivial way to satisfy conditions 1 \Gamma 4 of
Definition 2.4 is to assign "0" to all the row (column) nodes and then use positive colors on all
the column (row) nodes to satisfy condition 3. This ordering supports the tenet that use of bi-
partition/bi-coloring is never worse than one-sided calculation and that a substitution approach is
never worse than a direct approach (in principle).
Bi-coloring
The two combinatorial problems we face, corresponding to direct determination and determination
by substitution, can both be approached in the following way. First, permute and partition the
structure of J : ~
indicated in Figure 3.1. The construction of this
partition is crucial; however, we postpone that discussion until after we illustrate its' utility. Assume
I and
R
R

Figure

3.1: Possible partitions of the matrix ~
Second, define appropriate intersection graphs G I
R based on the partition [J C jJ R ]; a coloring
of G I
C yields a partition of a subset of the columns, GC , which defines matrix V . Matrix W is
defined by a partition of a subset of rows, GR , which is given by a coloring of G I
R . We call this
double coloring approach bi-coloring. The difference between the direct and substitution cases is
in how the intersection graphs, G I
R , are defined, and how the nonzeroes of J are extracted from
the pair (W T J; JV ).
3.1 Direct determination
In the direct case the intersection graph G I
C is defined: G I
9k such that J kr 6= 0; J ks 6= 0 and either (k; r) 2 JC or
The key point in the construction of graph G I
C , and why G I
C is distinguished from the usual
column intersection graph, is that columns r and s are said to intersect if and only if their nonzero
locations overlap, in part, in JC : i.e., columns r and s intersect if J kr \Delta J ks = 0 and either
The "transpose" of the procedure above is used to define G I
R ). Specifically, G I
I
R if nnz(row
If M is a matrix or a vector then "nnz(M)" is the number of nonzeroes in M
In this case the reason graph G I
R is distinguished from the usual row intersection graph is that
rows r and s are said to intersect if and only if their nonzero locations overlap, in part, in JR : rows
r and s intersect if J rk \Delta J
The bi-partition (GR ; GC ), induced by coloring of graphs G I
R and G I
C , is consistent with direct
determination of J . To see this consider a nonzero element (i; is in a group
of GC (corresponding to a color) with the property that no other column in GC has a nonzero in
row i: hence, element (i; j) can be directly determined. Analagously, consider a nonzero element
row r will be in a group of GR (corresponding to a color) with the property that no
other row in GR has a nonzero in column s: hence, element (r; s) can be directly determined. Since
every nonzero of J is covered, the result follows.
Example: Consider the example Jacobian matrix structure shown in Figure 3.2 with the partition
42 44

Figure

3.2: Example Partition
The graphs GC and GR formed by the algorithm outlined above are given in Figure 3.3. Coloring
GC requires 3 colors, while GR can be colored in two. Boolean matrices V and W can be formed in
the usual way: each column corresponds to a group (or color) and unit entries indicate column (or
row) membership in that group:
\Theta J
\Theta J 42 J 23 J 44 0
Clearly, all nonzero entries of J can be identified in either JV or W T J .
c
c 21212
r
r
r

Figure

3.3: Graphs GC and GR (direct approach)
3.2 Determination by substitution
The basic advantage of determination by substitution in conjunction with partition
is that sparser intersection graphs G I
R can be used. Sparser intersection graphs mean thinner
matrices V ,W which, in turn, result in reduced cost.
In the substitution case the intersection graph G I
C is defined: G I
9k such that J kr 6= 0; J ks 6= 0 and both (k; r) 2 JC , (k; s) 2 JC .
Note that the intersection graph G I
captures the notion of two columns intersecting if
there is overlap in nonzero structure in JC : columns r and s intersect if J kr \Delta J ks = 0 and both
some k. It is easy to see that E I
C is a subset of the set of edges used in
the direct determination case.
The "transpose" of the procedure above is used to define G I
R ). Specifically, G I
I
R if row i 2 JR and nnz(row
The intersection graph G I
R ) captures the notion of two rows intersecting if there is overlap
in nonzero structure in JR : rows r and s intersect if J rk \Delta J
for some k. It is easy to see that E I
R is a subset of the set of edges used in direct determination.
All the elements of J can be determined from (W T J; JV ) by a substitution process. This is
evident from the illustrations in Figure 3.4.

Figure

3.4 illustrates two of four possible nontrival types of partitions. In both cases it is clear
that nonzero elements in the section labelled "1" can be solved for directly - by the construction
process they will be in different groups. Nonzero elements in "2" can either be determined directly,
or will depend on elements in section "1". But elements in section "1" are already determined
(directly) and so, by substitution, elements in "2" can be determined after "1". Elements in
section "3" can then be determined, depending only on elements in "1" and "2", and so on until
the entire matrix is resolved.
JR
Figure

3.4: Substitution Orderings
Example. Consider again the example Jacobian matrix structure shown in Figure 3.2.
Column and row intersection graphs corresponding to substitution are given in Figure 3.5. Note
that GC is disconnected and requires 2 colors; GR is a simple chain and also requires 2 colors.
c
c 21112
r
r
r

Figure

3.5: Graphs GC and GR for substitution process
The coloring of GC and GR leads to the following matrices V , W and the resulting computation
of JV , W
\Theta 0
\Theta J
\Theta J 42 J 23 J 44 0
It is now easy to verify that all nonzeroes of J can be determined via substitution.
3.3 How to partition J.
We now consider the problem of obtaining a useful partition [J C jJ R ], and corresponding permutation
matrices P ,Q, as illustrated in Figure 3.1. A simple heuristic is proposed based on the
knowledge that the subsequent step, in both the direct and the substitution method, is to color
intersection graphs based on this partition.
Algorithm MNCO builds partition JC from bottom up, and partition JR from right to left.
At the k th major iteration either a new row is added to JC or a new column is added to JR : the
choice depends on considering a lower bound effect:
ae(J T
r
where ae(A) is the maximum number of nonzeroes in any row of matrix A, r is a row under consideration
to be added to JC , c is a column under consideration to be added to JR . Hence, the
number of colors needed to color G I
C is bounded below by ae(J C ); the number of colors needed to
color G I
R is bounded below by ae(J T
In algorithm MNCO, matrix C) is the submatrix of J defined by row indices R and
column indices C: M consists of rows and columns of J not yet assigned to either JC or JR .
Minimum Nonzero Count Ordering (MNCO)
1. Initialize
2. Find r 2 R with fewest nonzeros in M
3. Find c 2 C with fewest nonzeros in M
4. Repeat Until
if ae(J T
R=R-frg
else
C=C-fcg
repeat
Note that, upon completion, JR ; JC have been defined; the requisite permutation matrices are
implicitly defined by the ordering chosen in MNCO.
Bi-coloring performance
In this section we present results of numerical experiments. The work required to compute the
sparse Jacobian matrix is the work needed to compute (W T J; JV ) which, in turn, is proportional
to the work to evaluate the function F times the sum of the column dimensions of the boolean
matrices . The column dimension sum, is equal to the number
of colors used in the bi-coloring. In our experiments we compare the computed coloring numbers
required for the direct and substitution approaches. We also compute the number of colors required
by one-sided schemes: a column partition alone corresponds to the construction of V based on
coloring the column intersection graph of J , a row partition alone corresponds to the construction
of W based on coloring the row intersection graph of J . The latter case leads to the application of
the reverse mode of AD (alone), whereas the former case leads to use of the forward mode.
Both the direct and substitution methods require colorings of their respective pairs of intersection
graphs, G I
R . Many efficient graph coloring heuristics are available: in our experiments we
use the incidence degree (ID) ordering [3, 8].
We use three sources of test matrices: a linear programming testbed with results reported
in

Table

1 and summarized in Table 2; the Harwell-Boeing sparse matrix collection, with results
reported in Tables 3,4; self-generated m-by-n "grid matrices" with results given in Tables 5, 6.
A grid matrix is constructed in the following way. First, approximately
n of the columns are
chosen, spaced uniformly. Each chosen column is randomly assigned DENS \Delta m nonzeroes. Second,
approximately p m of the rows are chosen, spaced uniformly. Each chosen row is randomly assigned
nonzeroes. We vary DENS as recorded in Table 5.
For each problem we cite the dimensions of the matrix A and the number of nonzeros(nnz).
The experimental results we report are the number of colors required by our bi-coloring approach,
both direct and substitution, and the number of colors required by one-sided schemes.
Bi-coloring One-sided
Name m n nnz Direct Substitution column row
standata
stair 356 620 4021 36 29 36 36
blend 74 114 522
vtp.base 198 347
agg 488 615 2862 19 13 43 19
agg2 516 758 4740 26 21 49 43
agg3 516 758 4756 27 21 52 43
bore3d 233 334 1448 28 24 73 28
israel 174 316 2443 61
adlittle 56 138

Table

1: LP Constraint Matrices (http://www.netlib.org/lp/data/)

Table

2: Totals for LP Collection
Bi-coloring 1-sided Coloring
Name M N NNZ Direct Substitution column row
cannes 256 256 256 2916
cannes 268 268 268 1675
cannes
cannes 634 634 634 7228 28 21 28 28
cannes 715 715 715 6665 22
cannes 1054 1054 1054 12196 31 23
cannes 1072 1072 1072 12444
chemimp/impcolc 137 137 411 6 4 8 9
chemimp/impcold 425 425 1339 6 5 11 11
chemimp/impcole 225 225 1308 21 14
chemwest/west0067 67 67 294 9 7 9 12
chemwest/west0381 381 381 2157 12 9 29 50
chemwest/west0497 497 497 1727 22 19 28 55

Table

3: The Harwell-Boeing collection (ftp from orion.cerfacs.fr)
4.1 Observations
First, we observe that the bi-coloring approach is often a significant win over one-sided determi-
nation. Occasionally, the improvement is spectacular, e.g., "cannes 715". Improvement in the
Harwell-Boeing problems are generally more significant than on the LP collection in the sense that
bi-coloring significantly outperforms both one-sided possibilities. This is partially due to the fact
that the matrices in the LP collection are rectangular whereas the matrices in the Harwell-Boeing
collection are square: calculation of the nonzeroes of J from W T J alone can be quite attractive
when J has relatively few rows. The grid collection displays the advantage of bi-coloring to great
effect - grid matrices are ideal bi-coloring candidates.
In general the advantage of substitution over direct determination is not as great as the difference
between bi-coloring and one-sided determination. Nevertheless, fewer colors are almost always
needed and for expensive functions F this can be important. For most problems the gain is about
20% though it can approach 50%, e.g., "watt2".

Table

4: Totals for Harwell-Boeing Collection
Bi-coloring 1-sided coloring
M N DENS Direct Substitution column row
100 100 0.52 20 20 84 74
100 100 0.64 20 20 95 93
100 100 1.00 20 20 100 100
100 400 0.53
100 400 0.64
100 400 1.00

Table

5: Grid Matrices
4.2 Interface with ADOL-C
We have interfaced our coloring and substitution routines with the ADOL-C software. The C++
package ADOL-C [14] facilitates the evaluation of first and higher order derivatives of vector func-
tion, defined by programs written in C or C++.
We compare the time needed on a sample problem with respect to five approaches:
ffl AD/bi-coloring (direct)
ffl AD/bi-coloring (substitution)
ffl AD/column coloring (forward mode)
ffl AD/row coloring (reverse mode)
ffl FD (Sparse finite differencing based on column coloring)
The test function F we use is a simple nonlinear function: define
be the index set of nonzeroes in row i of the Jacobian matrix and define F

Table

Totals for Grid Matrices
The Jacobian matrix (and thus the sparsity pattern) is a 10-by-3 block version of Figure 2.3, i.e.,
3. Problem dimensions as indicated in Figure 4.1, were used in the
experiments.
Our results, portrayed in Figure 4.1, suggest the following order of execution time requirement
by different techniques:
Note that FD requires more time than AD=column even though the same coloring is used for both.
This is because the work estimate t V \Delta !(F ) is actually an upper bound on the work required by
the forward mode where t V is the number of columns of V . This bound is often loose in practise
finite-differencing since the subroutine to evaluate F is actually called
times.
Problem Size
Time
(in
seconds)
Performance graph of different sparse approaches
AD Column
AD Row
FD
AD-BiColoring-Direct
AD-BiColoring-Substitution

Figure

4.1: A comparison of different sparse techniques
Another interesting observation is that the reverse mode calculation (AD/Row) is about twice
as expensive as the forward calculation (AD/Column). This is noteworthy because in this example,
based on the structure Figure 2.1, the column dimensions of V and W are equal. This suggests
that it may be practical to weigh the cost of the forward calclulation of Jv versus the calculation
of w T J , where w; v are vectors. We comment further on this aspect in x7.
5 Substitution and round-off
In general, the substitution approach requires fewer colors and therefore is more efficient 3 , in
principle, than direct determination. However, there is a possibility of increased round-off error due
to the substitution process. In fact an analogous issue arises in the sparse Hessian approximation
context [4, 7, 16] where, indeeed, there is considerable cause for concern. The purpose of this
section is to examine this question in the AD context. The bottom line here is that there is less to
worry about in this case. In the sparse Hessian approximation case significant error growth occurs
when the finite-difference step size varies over the set of finite-difference directions; however, in our
current setting there this is not an issue since the "step size" is equal to unity in all cases.
First we consider the number of substitutions required to determine any nonzero of J from
are chosen using our substitution stategy. There is good news:
similar to the sparse Hessian approximation situation [4, 7, 16], the number of dependencies, or
substitutions, to resolve a nonzero of J can be bounded above by 1b(m
Theorem 5.1 Let OE be a bipartite cyclic p-coloring of G b (J). Then, OE corresponds to a substitution
determination of J and each unknown in J is dependent on at most m+n \Gamma 2 unknowns. Moreover,
it is possible to order the calculations so that the maximum number of substitutions is less than or
equal to 1
First, edges (nonzeros) with one end assigned color "0" can be determined directly: by
the definition of bi-coloring there will be no conflict. Second, every pair of positive colors induces a
forest (i.e., a collection of trees in G b (J)); therefore, the edges (nonzeros) in the induced forest can
be resolved via substitution [4]. Hence, all edges (nonzeros) can be resolved either directly or by
a substitution process and the worst case corresponds to a tree with m yielding an
upper bound of m+n \Gamma 2 substitutions. However, it is easy to see that the substitution calculations
can be ordered to halve the worst-case bound yielding at most 1
substitutions.
Next we develop an expression to bound the error in the computed Jacobian. Except for the
elements that can be resolved directly, the nonzero elements of the Jacobian matrix can be solved
for by considering each subgraph induced by 2 positive colors (directions), one color corresponding
to a subset of rows, one color corresponding to a subset of columns. Let us look at the subgraph
G p;q induced by colors p (columns), q (rows). Let z
q J , R q be the set of rows
colored q, and C p be the set of columns colored p, where
Let
denote the quantities computed via AD. Note that the errors introduced here are only due to the
automatic differentiation process and are typically very small.
In the solution process an element J ij is determined:
z
3 Of course a substitution method does incur the extra cost of performing the substitution calculation. However,
this can be done very efficiently and the subsequent cost is usually negligible.
depending on whether a column equation (of form Jv) or a row equation (of form w T J) is used.
Here, N(r i ) denotes set of neighbours of row i in G p;q , and N(c j ) denotes set of neighbours of
column j in G p;q .
Assume that J actual denotes the actual Jacobian matrix; hence,
J actual
J actual
or J actual
J actual
Define an error matrix ij be the difference
z
depending on the way element J ij was computed.
We take into account the effect of rounding errors by letting -
ij to be equal to ffl ij plus the
contribution from rounding errors due to use of the equation that determines J ij . We can now
again, depending on the way J ij is calculated.
Moreover, we let ffl max be the constant :
Note that ffl max has no contribution from step sizes, unlike results for finite- differencing [4, 16].
Theorem 5.2 If J is obtained by our substitution process then
Proof: From equation (5.1),
or
Let us assume, without loss of generality, that equation (5.2) holds. This implies a bound
But the same decomposition can be applied recursively to each E ik , and using Theorem 5.1, the
result follows.
There are two positive aspects to Theorem 5.2. First, unlike the sparse finite-difference substitution
method for Hessian matrices [4, 7, 16], there is no dependence on a variable "step size": in AD
the "step size" is effectively uniformly equal to unity. Second, there is no cumulative dependence
on nnz(J) but rather just on the matrix dimensions, m+ n. However, there is one unsatisfactory
aspect of the bound in Theorem 5.2: the bound is expressed in terms of ffl max , but ffl max is not
known to be restricted in magnitude. A similar situation arises in the [4, 7, 16]. Nevertheless, as
illustrated in the example discussed below, ffl max is usually modest in practise.
We conclude this section with a small experiment where we inspect final accuracies of the
computed Jacobian matrices. The test function F is a simple nonlinear function as described in
x4.2.
In

Table

7 "FD1" is the sparse finite difference computation [8] using a fixed stepsize
"FD2" refers to the sparse finite difference computation [8] using a variable stepsize: ff is uniformly
varied in the range [
ffl]. The column labelled "Rel error" records kERRk 2 , where the
nonzeros of ERR correspond to the nonzeros of J : for
computed
J actual
The general trends we observe are the following. First, similar to the results reported in [1]
for forward-mode direct determination, the Jacobian matrices determined by our bi-coloring/AD
approach are significantly and uniformly more accurate than the finite-difference approximations.
This is true for both direct determination and the substitution approach. Second, the direct
approach is uniformly more accurate than the substitution method; however, the Jacobian matrices
determined via substitution are sufficiently accurate for most purposes, achieving at least 10 digits
of accuracy and usually more. Finally, on these problem there is relatively little difference in
accuracy between the fixed step method and the variable step method. However, as
illustrated in [6], it is easy to construct examples where the variable step approach produces
unacceptable accuracy.
Direct Substitution FD1 FD2
size Rel error Rel error Rel error Rel error

Table

7: Errors (sample nonlinear problem)
6 Concluding remarks
We have proposed an effective way to compute a sparse Jacobian matrix, J , using automatic
differentiation. Our proposal uses a new graph technique, bi-coloring, to divide the differentiation
work between the two modes of automatic differentiation, forward and reverse. The forward mode
computes the product JV for a given matrix V ; the reverse mode computes the product W T J
for a given matrix W . We have suggested ways to choose thin matrices V; W so that the work
to compute the pair (W T J; JV ) is modest and so that the nonzero elements of J can be readily
extracted.
Our numerical results strongly support the view that bi-coloring/AD is superior to one-sided
computations (both AD and FD) with respect to the order of work required. Of course AD
approaches offer additional advantages over FD schemes: significantly better accuracy, no need to
heuristically determine a step size rule, and the sparsity pattern need not be determined a priori
[14].
Implicit in our approach is the assumption that the cost to compute Jv by forward mode AD
is equal to the cost of computing w T J by reverse mode AD, where v; w are vectors. This is true
in the order of magnitude sense - both computations take time proportional to !(F ) - but the
respective constants may differ widely. It may be pragmatic to estimate "weights" w
respect to a given AD tool, reflecting the relative costs of forward and reverse modes. It is very
easy to introduce weights into algorithm MNCO (x4:3) to heuristically solve a "weighted" problem,
is the number of row groups (or colors assigned to the rows), and - 2
is the number of column groups (or colors assigned to the columns). The heuristic MNCO can be
changed to address this problem by simply changing the conditional (LB) to:
R
Different weights produce different allocations of work between forward and reverse modes, skewed
to reflect the relative costs. For example, consider a 50-by-50 grid matrix with
x5), and let us vary the relative weighting of forward versus reverse mode: w
. The results of our weighted bi-coloring approach are given in Table 8.

Table

8: Weighted problem results
Finally, we note that the bi-coloring ideas can sometimes be used to efficiently determine relatively
dense Jacobian matrices provided structural information is known about the function F .
For example, suppose is a partially separable function,
t, and each component function F i typically depends on only a few
components of x. Clearly each Jacobian function J i is sparse while the summation,
may or may not be sparse depending on the sparsity patterns. However, if we define a "stacked"
function ~
F ,
~
then the Jacobian of ~
F is
~
J is sparse and the bi-coloring/AD technique can be applied to ~
J, possibly yielding a
dramatic decrease in cost. Specifically, if J is dense (a possibility) then the work to compute J
without exploiting structure is n \Delta !(F ) whereas the cost of computing ~
J via bi-coloring/AD is
approximately
J)) is the minimum number of colors required
for a bipartite cyclic coloring of graph G b ( ~
J). Typically, - c (G b ( ~
n. The idea of applying the
bi-coloring/AD technique in a structured way is not restricted to partially separable functions [9].

Acknowledgements

We are very grateful to Andreas Griewank, his student Jean Utke, and his colleague David Juedes
for helping us with the use of ADOL-C.



--R



New methods to color the vertices of a graph
The cyclic coloring problem and estimation of sparse Hessian matrices




Structure and efficient Jacobian calculation
On the estimation of sparse Jacobian matrices
Direct calculation of Newton steps without accumulating Jacobians



Computing a sparse Jacobian matrix by rows and columns
On the estimation of sparse Hessian matrices
--TR

--CTR
Shahadat Hossain , Trond Steihaug, Sparsity issues in the computation of Jacobian matrices, Proceedings of the 2002 international symposium on Symbolic and algebraic computation, p.123-130, July 07-10, 2002, Lille, France
Dominique Villard , Michael B. Monagan, ADrien: an implementation of automatic differentiation in Maple, Proceedings of the 1999 international symposium on Symbolic and algebraic computation, p.221-228, July 28-31, 1999, Vancouver, British Columbia, Canada
L. F. Shampine , Robert Ketzscher , Shaun A. Forth, Using AD to solve BVPs in MATLAB, ACM Transactions on Mathematical Software (TOMS), v.31 n.1, p.79-94, March 2005
Shaun A. Forth, An efficient overloaded implementation of forward mode automatic differentiation in MATLAB, ACM Transactions on Mathematical Software (TOMS), v.32 n.2, p.195-222, June 2006
Thomas F. Coleman , Arun Verma, ADMIT-1: automatic differentiation and MATLAB interface toolbox, ACM Transactions on Mathematical Software (TOMS), v.26 n.1, p.150-175, March 2000
sparse finite differencing;NP-complete problems;ADOL-C;automatic differentiation;nonlinear systems of equations;bicoloring;computational differentiation;sparse Jacobian matrices;partition problem;nonlinear least squares;graph coloring