% Template article for preprint document class `elsart'
% SP 2001/01/05

% Modified 2007/01/09 for CTW 2010
% Please find the original file at http://www.elsevier.com/

\documentclass{ctwart}
\linespread{0.80}
\usepackage{amssymb}

\begin{document}

\begin{frontmatter}

\title{Graph Models and their Efficient Implementation for Sparse Jacobian 
Matrix Determination\thanksref{ac}}
\thanks[ac]{This research was supported in part by the Natural Sciences and Engineering Research Council of Canada (NSERC) and the Research Council of Norway (NFR).}

\author[sh]{Shahadat Hossain\corauthref{ca}} and \author[ts]{Trond Steihaug}

\address[sh]{Department of Mathematics and Computer Science, University of Lethbridge, Canada}
\address[ts]{Department of Informatics, University of Bergen, Norway}
 \corauth[ca]{Corresponding author.}
 \begin{abstract}
Algorithms for solving  large-scale combinatorial scientific computing problems arising in sparse or 
otherwise structured  
matrix computation  are often represented by appropriate graph 
models and sometimes the same problem can be 
formulated in more than one graph models with similar asymptotic computational complexity. The relative merits of 
different graph models for the same problem can then be expressed in terms of factors such as 
generality of the model and ease of computer implementation.  In this note we briefly review the contemporary graph formulations 
for large-scale sparse Jacobian matrix determination problem (JMDP) and 
suggest the pattern graph model which can be viewed as a unifying framework 
for the unidirectional and bidirectional approaches for JMDP.
Due to 
the irregular memory access pattern combined with low floating point calculations 
relative to the volume of data movements  the actual 
running time of sparse matrix and graph algorithms may achieve only a small fraction 
of the theoretical performance.  We proffer the use of
array-based  data structures as  the basic  
building-blocks for efficient implementation of fundamental 
graph algorithms  on modern cache-based computer architectures. Numerical 
results comparing our implementation (DSJM toolkit) with ColPack \cite{Pothen2009a}
is given.
  \end{abstract}

\begin{keyword}
Sparse Matrix Data Structures \sep Intersection Graph \sep Bipartite Graph
\sep Hypergraph \sep Pattern Graph 
\end{keyword}

\end{frontmatter}

% main text
\section{Introduction}
\label{intro}
We consider the problem of determining the Jacobian matrix $F^\prime(x)$ of a mapping 
$F: \Re^n \mapsto \Re^m$.
 In this paper graphs are undirected.   
The colon notation of
  \cite{Golub96a} is used to  specify sections of a matrix. 
  The $(i,j)$th entry is denoted by $A(i,j)$ or  $a_{ij}$. When a matrix is displayed the nonzero entries are explicitly shown while a blank or 
  a '$0$' marks a zero entry. We assume that the sparsity pattern of matrix $F^\prime(x)$
  is known a priori and that it is computationally more economical to compute the entire 
  column of  $F^\prime(x)$  than computing individual entries. Using  differences 
  the $j$th column of the Jacobian matrix may be approximated as 
  \begin{equation}
\left .\frac{\partial F(x+ts)}{\partial t}\right|_{t=0}  = F^\prime(x)s
\approx As=\frac{1}{\varepsilon}[F(x+ \varepsilon s) - {F(x)}]\equiv b
\label{fd}
\end{equation}
 with one extra function evaluation. Also 
Algorithmic Differentiation (AD) \cite{Griewank2008a} forward mode gives $b=F^\prime(x)s$ at a cost 
which is a small multiple of the cost of one function evaluation. The Jacobian matrix 
determination problem can be stated as below.
\begin{prob}[JMDP]
Obtain vectors
$s_i \in {\Re}^n, i = 1, \ldots, p \mbox{ and }  w_j \in {\Re}^m, j = 1, \ldots, q\mbox{ with } p+q \mbox{ minimized}$ such that  the products
$(b_k =  As_k, \;  j=1,\cdots,p$ or  $B = AS)$  and  
$(c_k^T = w_k^T A, \;  k=1,\cdots,q$  or  $C^T = W^TA)$
determine the matrix  $A$ uniquely.
\label{jdp}
\end{prob}
In absence of any sparsity information, one may use the Cartesian basis vectors 
$e_i, i = 1, \ldots, n$ in (\ref{fd}) using $n$ extra function evaluations. However,
if the columns are structurally orthogonal i.e. no two columns have nonzero entries 
in the same row position only one extra function evaluation 
\[\begin{array}{c}F^\prime_j + F^\prime_k \approx A(:,j) +  A(:,k) =
\frac{1}{\varepsilon}[F(x+ \varepsilon (e_j + e _k)) - {F(x)} ] \end{array} \]
is sufficient to read off  the nonzero entries from the product $b = As$. If the 
columns can be partitioned into $p$ structurally orthogonal 
groups then the  Jacobian matrix is {\em directly determined} from the 
compressed matrix $B = AS$. Similarly,  the rows can be 
partitioned into $q$ structurally orthogonal groups and the 
Jacobian can be directly determined from $C^T = W^T A$ using the reverse mode 
of AD.  The problem of finding minimum cardinality orthogonal column (or row) 
partition of matrix $A$ can be solved as different vertex coloring problems (which in
general are NP-hard and therefore are solved by  heuristics)  of  suitable
graph(s) associated with $A$.    
\section{Graph Models and  Computer Implementation}
With regard to the choice of the graph model for the partitioning problem we expect  that
 the graph formulation
\begin{enumerate}
\item retains exploitable matrix structures,
\item enables  efficient implementation of pertinent algorithms, and 
\item generic enough to encapsulate the combined row-and-column 
determination as in Problem JMDP. 
\end{enumerate}
 
 The {\em intersection graph} of the columns of  $A$ is
denoted
 by $G(A) = (V,E)$ where corresponding to $A(:,j), j = 1, 2, \ldots, n$, there is a vertex $v_j \in V$
and $\{v_j,v_l\} \in E$ if and only if  $A(:,j)$ and $A(:,l), l
\neq j$ have nonzero elements in the same row position.   Then an orthogonal 
partition of the columns of $A$ is equivalent to a coloring $\phi$ of the vertices of 
$G(A)$ such that  $\phi(u) \ne \phi(v)$ if and only if $\{u,v\} \in E$ \cite{Coleman1983a}. Unfortunately, 
the column intersection graph of $A$ or $A^T$ is unable to expose all the exploitable
matrix sparsity. Alternative formulations define graphs based on column segments 
to allow for better utilization of available structure or sparsity \cite{HossainS08}.  However, these 
alternative formulations apply to either column direction or row direction but not 
to a combination of row and column directions -- henceforth {\em bidirectional determination}.     


 The {\em bipartite graph} associated with  matrix $A$ is
denoted
 by $G_b(A) = (V_c \cup V_r, E)$ where corresponding to $A(:,j), j = 1, 2, \ldots, n$, there is a 
 column vertex $v_j \in V_c$ and  corresponding to $A(i,:), i = 1, 2, \ldots, m$, there is a 
 row vertex $v_i \in V_r$
and $\{v_i,v_j\} \in E$ if and only if  $a_{ij} \ne 0, i = 1,2, \ldots, m,
 j = 1,2, \ldots, m$. The bipartite graph model has been proposed independently by Coleman and Verma  
 \cite{Coleman1998a} and Hossain and Steihaug \cite{Hossain1998a} in connection 
 with bidirectional determination of Jacobian matrices. In \cite{GMP2005a} the bipartite graph model
 has been considered for  column partitioning.  The zero-nonzero structure of the 
 underlying matrix is accurately represented  in  its bipartite graph making it 
 a natural logical data structure for bidirectional determination. On the other hand, for the unidirectional determination the model  is ``asymmetrical'' in the sense that  it contains extraneous information.  This difficulty  is manifested in \cite{GMP2005a} where the unidirectional determination posed 
 as distance-$2$ coloring needed the qualification ``partial'' as only one set of 
 vertices are colored. 
 
A {\em hypergraph} is a graph in which edges are generalized as {\em hyperedges} where a hyperedge may connect more than vertices. 
The hypergraph model presented here 
is more general than the ones considered in \cite{GMP2005a}. The  hypergraph 
$H(A) = (V, E)$ associated with the matrix $A$ has the vertex set
\[V =  \{v_i  | \exists k \mbox{ for which } a_{ik} \ne 0, i = 1,2, \ldots, m \} \cup 
   \{v_j  | \exists k \mbox{ for which } a_{kj} \ne 0, j = 1,2, \ldots, n,\}   
 \]
 and  corresponding to each row  $i$ and  each column $j$ the hyperedges $e_i \in E$ 
 and  $e_j \in E$, respectively,  are defined by  
 \[ e_i = \{ v_k |  a_{ik} \ne 0\} \mbox{ and }  e_j = \{ v_k |  a_{kj} \ne 0\}.\]
 
A {\em lateral neighbor} of $a_{ij} \ne 0$ is a nonzero $a_{ij^\prime} \ne 0$ in row $i$ 
of $A$ such that $j^\prime - j$ is the smallest if $j^\prime > j$ or such that 
 $j - j^\prime$ is the smallest if $j > j^\prime$ among all  such indices $j^\prime$ in row $i$. 
Vertical neighbors can be interpreted  in an analogous way with the roles of $i$ and $j$ interchanged.
 The {\em sparsity-pattern graph} 
(or simply the pattern graph)  associated with $A$ is  $G_{\mathcal{P}}(A) = (V, E)$, 
where \[V = \{v_{ij}:a_{ij} \neq 0, i = 1,2, \ldots, m, j = 1,2,\ldots, n\}\] and \[ \{v_{ij},v_{i^\prime j^\prime}\} \in E
\mbox{ if } a_{ij} \mbox{ and  } a_{i^\prime j^\prime}  \mbox{ are lateral or vertical neighbors.}\]    

 
 An unknown $a_{ij}$ is said to be {\em covered} (by $S$ or $W$) if it can be uniquely solved in 
\[ \hat{S_i}^T A(i,\mathcal{J}_i)^T = B(i,:)^T \mbox{ or }  A(\mathcal{I}_j, j)  \hat{W_j} = C(:,j). \]
where $\hat{S_i} (\hat{W_j})$ is the submatrix of $S (W)$ corresponding to the nonzero 
entries in row $i$ (column $j$) indicated by $\mathcal{J}_i (\mathcal{I}_j)$. The matrices $S$ and $W$ are said to constitute a cover for $A$ if each $a_{ij} \ne 0$ is 
covered. A cover is a {\em direct cover} if $A$ can be  determined  directly from the cover.
Given a  mapping   $\Phi: V \mapsto \mathcal{S} \cup \mathcal{W}$ where 
 $\mathcal{S} = \{\mathcal{S}_1, \mathcal{S}_2, \ldots, \mathcal{S}_p\}, \mathcal{W} = \{\mathcal{W}_1, \mathcal{W}_2, \ldots, \mathcal{W}_q \}$
 we  define  matrices  $S \in \{0,1\}^{n \times p}$ and $W \in \{0,1\}^{m \times q},$
\[ S(:,k) = \sum_j e_j, \mathcal{S}_k = \Phi(v_{ij}) \mbox{ and } W(:,l) = \sum_i e_i,  \mathcal{W}_l = \Phi(i).\]
The mapping 
$\Phi: V \mapsto \mathcal{S} \cup \mathcal{W}$ is  said to yield a cover for $A$ if  the matrices $S$ and $W$  
constitute a cover for $A$. Denote by $ u_{ij} \stackrel{\ge 1}{\sim} u_{i^\prime  j^\prime}$
a path of length at least $1$.
\begin{thm}
Let $G_{\mathcal{P}}(A) = (V, E)$ be the pattern graph associated with $A$. Define the 
 mapping 
 $\Phi: V \mapsto \mathcal{S} \cup \mathcal{W}$ where 
 $\mathcal{S} = \{\mathcal{S}_1, \mathcal{S}_2, \ldots, \mathcal{S}_p\}, 
\mathcal{W} = \{\mathcal{W}_1, \mathcal{W}_2, \ldots, \mathcal{W}_q \}$ such that for each $v_{ij} \in V,$
  \begin{enumerate}
 \item[] EITHER
\item \begin{enumerate}
 \item   $ v_{ij} \stackrel{\ge 1}{\sim} v_{i  j^\prime},   j \ne j^\prime 
  \mbox{ implies }k \ne k^\prime \mbox{ where } \Phi(v_{ij}) =  \mathcal{S}_k,  \Phi(v_{ij^\prime}) 
  = \mathcal{S}_{k^\prime} \mbox{ and }$
\item  $v_{ij} \stackrel{\ge 1}{\sim} v_{i  j^\prime} \stackrel{\ge 1}{\sim} v_{i^\prime  j^\prime},  i \ne i^\prime, j \ne j^\prime 
 \mbox{ implies } k \ne k^\prime \mbox{ where } \Phi(v_{ij}) =  \mathcal{S}_k,  \Phi(v_{i^\prime j^\prime}) 
  = \mathcal{S}_{k^\prime} $
    \end{enumerate}
    \item[]  OR
  \item \begin{enumerate}
 \item   $ v_{ij} \stackrel{\ge 1}{\sim} v_{i ^\prime j},   i \ne i^\prime 
  \mbox{ implies }k \ne k^\prime \mbox{ where } \Phi(v_{ij}) =  \mathcal{W}_k,  \Phi(v_{i^\prime j}) 
  = \mathcal{W}_{k^\prime} \mbox{ and }$
\item  $v_{ij} \stackrel{\ge 1}{\sim} v_{i^\prime  j} \stackrel{\ge 1}{\sim} v_{i^\prime  j^\prime},  i \ne i^\prime, j \ne j^\prime 
 \mbox{ implies } k \ne k^\prime \mbox{ where } \Phi(v_{ij}) =  \mathcal{W}_k,  \Phi(v_{i^\prime j^\prime}) 
  = \mathcal{W}_{k^\prime}. $
    \end{enumerate}
    \end{enumerate}
    Then matrices $S$ and $W$ constitute a direct cover for $A$.
    \end{thm}
One of the strengths of the pattern graph model is that the  result given above can be specialized to unidirectional 
and column segments determinations without changing the graph.

The sparsity pattern of matrix $A$ can be efficiently represented using  two arrays: array \verb"colind" that stores the column indices of the nonzero entries row-by-row, and array \verb"rowptr" that
 contains the index of the first nonzero element of each row of the sparse matrix stored in
 \verb"colind" array. For easy access to the adjacent vertices sparsity pattern of 
 $A^T$ is explicitly stored using analogous arrays \verb"rowind"  and \verb"colptr" . This storage is of  order 
 $\Theta(\max(n,m,nnz))$ which meets the design strategies  for sparse linear algebra implementation
 \cite{GilbertRS06}. Next 
 we consider the computations on this data structure relevant to algorithms for Problem JMDP.  In many graph coloring heuristics on static graphs a frequently executed 
 operation is to find the neighbors of a given vertex $v_j$. Assuming $v_j$ a column vertex  this information is obtained easily as 
 \[\left\{\{v_j, v_i\}| i = \mbox{  rowind(k), k = colptr(j) : colptr(j+1)-1 }  \right\}.\]
 This computation can be performed independent of  the graph models discussed above. Further,
 the array representation ensures better cache performance and  the computational 
 cost of the operation is proportional to the size of the data accessed and the 
 number of nonzero arithmetic operations. Clearly, it is not necessary 
 to compute the intersection graphs explicitly as advocated in \cite{GMP2005a}. Indeed,
 almost all the well-known coloring heuristics for  JMDP can be implemented in 
 time proportional to $\sum_{i=1}^{m} \rho_i^2$ where $\rho_i$ denotes 
 the number of nonzero entries in row $i$ of $A$. As an indication of the 
 efficiency of the data structure proposed here we report in the table below  
 the  timing experiments 
 for incidence degree order (IDO) coloring where 
 \begin{table}[hbt]
\begin{tabular}{|crrr|rr|rr|}
\hline
{\em Matrix} &$m$ &$n$ &$nnz$ &\multicolumn{2}{|c|}{ColPack}\cite{Pothen2009a}& \multicolumn{2}{|c|}{DSJM}\\ \hline
&&&&ot&ct&ot&ct\\ \cline{5-8} 
lpcreb &   9648 &  77137 &  260785 & 13.9 & 1.49 & 2.46 & 1  \\
 lpcred &   8926 &  73948 &  246614 & 14.79 & 1.48 & 2.46 & 1.01 \\
  lpfit2d &     25 &  10524 &  129042 & 64.19 & 24.94 & 16.32 & 17.2 \\
   lpken18 & 105127 & 154699 &  358171 & 15.92 & 0.63 & 1.47 & 0.48\\
   lposa07 &1118 &25067 &144812  & 161.76 & 28.34 & 40.84 & 18.35 \\ \hline
\end{tabular}
 \label{t1}
\end{table}
\verb"ot"  and \verb"ct" represent ordering and coloring times, respectively. 
\verb"ColPack" implements IDO coloring with  partial distance-$2$ coloring scheme
on bipartite graph while \verb"DSJM" implements IDO coloring using the data structure described 
here using column intersection graph. The times reported here do not include 
data structure set up times. The sparse matrices are obtained from University of 
Florida Sparse Matrix Collection.  

\begin{thebibliography}{00}

\bibitem{Coleman1983a}
T.~F. Coleman and J.~J. Mor{\'e}.
\newblock Estimation of sparse {J}acobian matrices and graph coloring problems.
\newblock {\em SIAM J. Numer. Anal.}, 20(1):187--209, 1983.

\bibitem{Coleman1998a}
T.~F. Coleman and A.~Verma.
\newblock The efficient computation of sparse {J}acobian matrices using
  automatic differentiation.
\newblock {\em SIAM J. Sci. Comput.}, 19(4):1210--1233, 1998.

\bibitem{GMP2005a}
A.~H. Gebremedhin, F.~Manne, and A.~Pothen.
\newblock {W}hat {C}olor {I}s {Y}our {J}acobian? {G}raph {C}oloring {F}or
  {C}omputing {D}erivatives.
\newblock {\em {SIAM} {R}eview}, 47(4):629--705.

\bibitem{Pothen2009a}
A.~H. Gebremedhin, A.~Tarafdar, D.~Nguyen, and A.~Pothen.
\newblock   ColPack.
\newblock http://www.cs.odu.edu/\verb"~"dnguyen/dox/colpack/html/ (accessed May 2009)

\bibitem{GilbertRS06}
J.~R. Gilbert, S.~Reinhardt, and V.~Shah.
\newblock High-performance graph algorithms from parallel sparse matrices.
\newblock In B.~K{\aa}gstr{\"o}m, E.~Elmroth, J.~Dongarra, and J.~Wasniewski,
  editors, {\em PARA}, volume 4699 of {\em Lecture Notes in Computer Science},
  pages 260--269. Springer, 2006.

\bibitem{Golub96a}
G.~H. {Golub} and C.~F. {Van Loan}.
\newblock {\em Matrix Computations}.
\newblock Johns Hopkins University Press, 3rd edition, 1996.

\bibitem{Griewank2008a}
A.~Griewank, Andrea Walther.
\newblock {\em Evaluating Derivatives: Principles and Techniques of Algorithmic
  Differentiation}.
\newblock Number~19 in Frontiers in Appl. Math. SIAM, Philadelphia, Penn.,
  2008.

\bibitem{HossainS08}
 Shahadat Hossain and Trond Steihaug.         
  \newblock Graph coloring in the estimation of sparse derivative matrices:
               Instances and applications,
  \newblock{\em Discrete Applied Mathematics}, 156(2):280--288, 2008.
  
\bibitem{Hossain1998a}
A.~S. Hossain and T.~Steihaug.
\newblock Computing a sparse {J}acobian matrix by rows and columns.
\newblock {\em Optimization Methods and Software}, 10:33--48, 1998.

\end{thebibliography}

\end{document}
