Triangle counting via sparse linear algebra¶

This example considers different ways to count triangles in a graph. It illustrates some of the constructs available in Intrepydd, although for this particular problem, a pure Numpy/Scipy version is likely faster in the sequential case. (On a multicore system, an Intrepydd implementation may be faster or more scalable.)

Background: Counting triangles in a social network¶

A social network may be modeled as an undirected graph, like the one shown below.

The nodes (or vertices) of this graph, shown as numbered circles, might represent people, and the edges (or links connecting them) might represent who is friends with whom. In this case, person 0 is friends with all the "odd birds" of this network, persons 1, 3, and 5, but has no direct connection to persons 2 and 4.

One type of analysis one might perform on such a graph is counting triangles, that is, the number of relationships of the form $a$ knows $b$, $b$ knows $c$, and $c$ knows $a$. In the graph shown above, there are two such triangles: (0, 1, 3) and (0, 3, 5).

A conventional graph-based approach. A commonly used Python-based library for graph computations is NetworkX. While it turns out to be much slower than the method described in this notebook, if you are interested in seeing how to use NetworkX, see our supplemental notebook.

A matrix-based method for counting triangles¶

Adjacency matrix. Let $A$ be the adjacency matrix representation of the graph, defined as follows. The entries of $A$ are either 0 or 1; and $a_{i,j}$ equals 1 if and only if there is an edge connecting nodes $i$ and $j$. For instance, for the graph shown above,

$$A = \begin{bmatrix} 0 & 1 & 0 & 1 & 0 & 1 \\ 1 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 1 & 1 & 1 & 0 & 1 & 1 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 1 & 0 & 0 \end{bmatrix}.$$

Observe that the relationships are symmetric. For instance, 0 and 1 are mutually connected; therefore, both $a_{0,1}$ and $a_{1, 0}$ equal 1, and in general, $A = A^T$.

Matrix-based graph analysis via linear algebra. Given such a matrix-based representation of a graph, here is one way to count triangles using linear algebra operations.

First, let $A \cdot B$ denote matrix multiplication. That is, $C = A \cdot B$ means $c_{i,j} = \sum_k a_{i,k} b_{k, j}$.

Next, let $A \odot B$ denote elementwise multiplication. That is, $E = A \odot B$ means $e_{i, j} = a_{i, j} b_{i, j}$.

Then, here is a two-step method to compute the number of triangles, which we'll denote as $t(A)$:

$$\begin{eqnarray} C & = & (A \cdot A) \odot A \\ t(A) & = & \frac{1}{6} \sum_{i, j} c_{i,j}. \end{eqnarray}$$

The first step computes a "count" matrix $C$. Each element, $c_{i,j}$, counts the number of triangles in which both $i$ and $j$ appear. You can see this fact algebraically; the scalar formula for $c_{i,j}$ is given by

$$c_{i,j} = \underbrace{\sum_k a_{i,k} \cdot a_{k, j}}_{A \cdot A} \cdot \underbrace{a_{i,j}}_{\odot A}$$

Observe that the summand involves the product $a_{i,k} a_{k, j} a_{i,j}$, which is nonzero only if all of the $(i,k)$, $(k,j)$, and $(i,j)$ edges exist.

For the example shown above, $c_{0, 1} = c_{1,0} = 1$ since there is only one triangle that uses the edge $(0, 1)$, whereas $c_{0, 3} = c_{3, 0} = 2$ because the edge $(0, 3)$ appears in two triangles.

The second step sums all the elements of $C$. However, the sum alone will overcount the number of unique triangles by a factor of six (6), hence the additional factor of $\frac{1}{6}$. (Why?)

Instead of summing all the entries of $A$, one can exploit symmetry and consider just the upper- or lower-triangle, which is a good additional optimization you might try once you've finished the main exercises contained herein.

The function, count_triangles(A, ret_C), below, implements the above formula. That is, given a symmetric Numpy array A representing the adjacency matrix of a graph, this function will return the number of triangles. It assumes the input matrix A is a (dense!) Numpy array.

The ret_C flag is a boolean, and when True, the function returns both the triangle count and the matrix $C$; otherwise, it returns just the triangle count.

The following cell does some extensive testing based on random sparse graphs.

The IMDB actor network dataset¶

Later, we will apply the matrix-based triangle counting method to a real "social" network, namely, the graph of actors who have appeared in the same movie. That is, let each actor be a node in the graph, with an edge between two actor nodes $i$ and $j$ if they appeared in the same movie together. If there are $n$ actors, then this graph can be represented by an $n \times n$ adjacency matrix, where $a_{i,j}=1$ if and only if $i$ and $j$ appeared in at least one movie together, and 0 otherwise. (We will only allow 0 or 1 values in the matrix, effectively ignoring the number of movies in which the actors "co-occurred.)

The dataset in this problem uses data collected on a crawl of the Top 250 movies on the Internet Movie Database (circa 2012). If you are interested in the raw data to construct the sparse matrix corresponding to the graph, refer to the NetworkX-based supplemental notebook.

Observation: sparsity. The matrix of actor co-occurrences is quite sparse. There are $12,975$ actors, which means as many as $\approx 12,\!975^2 / 2 \approx 84$ million possible co-occurences; yet, the matrix (A_casts) only appears to have $1.2$ million nonzeros, which is about $600,\!000$ observed co-occurrences, or less than 1% of the possible co-occurrences.

A sparse matrix-based approach¶

The first matrix-based example we gave above stores the graph as a dense matrix. However, as noted above, the graph itself is actually quite sparse. Thus, any method that scales like $\mathcal{O}(n^2)$ or worse, where $n$ is the number of actors, will be much slower than a method that exploits sparsity.

Luckily, Numpy/Scipy has built-in support for sparse matrices and basic operations on them. Therefore, we can use the simple implementation below to implement the counting formula, and trust that Numpy/Scipy works reasonably well if the input matrix is sparse.

Observation. In experiments we have done on conventional desktop machines, this sparse matrix-based implementation is 10-30x faster than a NetworkX-based implementation.

An Intrepydd version: Domain-specific wrappers¶

The linear algebraic algorithm has a "straightforward" translation to Intrepydd's domain-specific wrappers. In this example, the calc_C() Intrepydd kernel defined below implements the Numpy/Scipy statement, C = (A @ A).multiply(A), where A and C are sparse matrices. In the current version of Intrepydd, sparse matrix operations are wrappers for the CombBLAS library.

Take a moment to read this implementation, and observe the following.

• In the current version of Intrepydd, it is not possible to pass full sparse matrix objects between Python and Intrepydd. One workaround is to "extract" the underlying arrays representing the sparse matrix and pass those instead. Here, we assume the input A and the output C are compressed sparse row (CSR) matrices, which are represented by an array of floating-point values, integer column indices, and an integer array of row pointers.
• Sparse matrices in Intrepydd, at present, must store the values as floating-point. By contrast, the array A_casts constructed above has only integer values. That will require the caller to do an extra conversion, which might affect execution time.
• The Intrepydd call csr_to_spm() translates from these CSR arrays to an Intrepydd sparse matrix object. Refer to the Intrepydd Guide function list for additional details.
• Given an Intrepydd sparse matrix object A, the Intrepydd call spmm(A, A) corresponds to A @ A in Numpy/Scipy.
• The Intrepydd call spm_mul(A, B) corresponds to A.multiply(B).
• The Intrepydd call spm_to_csr() translates from an Intrepydd sparse matrix object back to an array representation for the caller. These are passed-by-reference arguments, which the caller of calc_C must supply. (See the modified Python function, count_triangles_spmat__pydd(), a couple code cells below.) Furthermore, this operation might also affect overall execution time, just as csr_to_spm() can.

The above defined the Intrepydd kernel; here is how it may be called from a Python driver function (count_triangles_spmat__pydd). Pay particular attention to how one gets the array representation for a Numpy/Scipy CSR object, and how one can construct the output.

Regarding the output, this implementation exploits the fact that for C = A.multiply(A @ A), C must have the same sparsity pattern as A. (Can you see why?) Therefore, when the driver allocates space for C (c_vals, c_inds, and c_ptrs), it is sufficient for these to have the same sizes as those for A.

Observation. The time is a little bit slower than the pure Numpy/Scipy version. Again, as noted previously, there are extra overheads in the current version of Intrepydd for input/output data transfers due to csr_to_spm() and spm_to_csr() calls.

Another Intrepydd version: "lowering" the implementation¶

The current version of Intrepydd is translated into C++ and then compiled. So, you can also try to write "low-level" implementations that look like what you might have written in C++, and hope to get C++-level performance from it. The following implementation is an example of "lowering" the original implementation.

In this case, lowering makes it possible to optimize across function calls. That is, recall the linear algebraic operation that triangle counting expresses, written here in a slightly different and more suggestive way (since $A \odot B = B \odot A$):

$$c_{i,j} = \underbrace{a_{i,j} \, \cdot}_{A \odot} \underbrace{\sum_k a_{i,k} \cdot a_{k, j}}_{A \cdot A}.$$

In this form, it should be clear that $C$ only has a nonzero output $c_{i,j}$ when $A$ has a nonzero element $a_{i,j}$. Therefore, rather than compute the entire product $B = A \cdot A$ first, we can compute only those $b_{i,j}$ where $a_{i,j} \neq 0$. The computation of a single $b_{i,j}$ is a sparse dot product between a sparse row $i$ of $A$ and a sparse column $j$ of $A$. And since $A$ is symmetric ($A = A^T$), sparse rows and columns are "interchangeable."

Of course, it is not clear if such an implementation will be faster or not because there are many ways to implement $A \cdot A$ to save computation and data movement, but from the point of view of learning Intrepydd, it is instructive to try.

Aside. This formula also explains why one expects $C$ to have the same sparsity pattern as $A$, as noted previously.

In addition, lowering confers one more potential benefit in this case. Recall that Intrepydd sparse matrix objects must store their values in floating-point; by contrast, if we lower and do not use the sparse interface at all, we can change the implementation to use integer types, too.

Logical operations. At present, logical operations (or, and, not) are not supported in Intrepydd. Therefore, the sparse dot product code cannot use the more compact, has_nzs_left = (ki < ni) and (kj < nj).

Observation. This Intrepydd version is not faster than the earlier Intrepydd version or the native Numpy/Scipy versions. However, having a lowered version means you can try even more optimizations, including parallelization! (All the other versions have been purely sequential.)

Fin! We hope this notebook has been a useful tutorial of how to use Intrepydd's sparse matrix operations, in the context of implementing graph algorithms that exploits the duality between a sparse matrix and a sparse graph. In particular, if you are working on a graph analysis problem, you might consider whether there is a matrix version and try to implement that using the APIs and techniques described above.