# Parallelize the construction of sparse matrices

I'm trying to parallelize a few constructions of chain complexes (given as a list of sparse matrices). I have a function that for any column (or row) computes the list of nonzero entries in that column (row), but several entries can appear at the same place $$(i,j)$$ in the matrix (their values should be summed).

For example, given a simplicial complex as a list bases (simplices of given dimension), I write:

chCx[bases_] := Module[
{dim=Length@bases, dims=Length/@bases, basesk, baseskk, bdrs={}, bdr, x},
basesk =Association@Table[bases[[1,i]]->i,{i,dims[[1]]}];
Do[baseskk=Association@Table[bases[[k,i]]->i,{i,dims[[k]]}];
...
AppendTo[bdrs,bdr]; basesk=baseskk;, {k,2,dim}]; bdrs];


Then I replace ... with 4 different commands, to get functions chCx, chCx0, chCx1, chCx2:

bdr = SparseArray[{},Length/@{basesk,baseskk}];
Do[bdr[[basesk[Delete[s,{{r}}]],baseskk[s]]]+=(-1)^(r+1),{r,1,k},{s,Keys@baseskk}];


and

bdr=SparseArray[{},Length/@{basesk,baseskk}]; SetSharedVariable[bdr]; DistributeDefinitions[basesk,baseskk];
ParallelDo[bdr[[basesk[Delete[s,{{r}}]],baseskk[s]]]+=(-1)^(r+1), {r,1,k},{s,Keys@baseskk}];


and

bdr = ParallelCombine[ Module[{bdri=SparseArray[{},Length/@{basesk,baseskk}]},
Do[bdri[[basesk[Delete[s,{{r}}]],baseskk[s]]]+=(-1)^(r+1), {r,1,k},{s,#}]; bdri]&, Keys@baseskk, Plus];


and

bdr = SparseArray[Flatten[ParallelTable[{basesk[Delete[s,{{r}}]],baseskk[s]}->(-1)^(r+1),
{r,1,k},{s,Keys@baseskk}],1],Length/@{basesk,baseskk}];


When I use these commands with

n=12; bases=Table[Subsets[Range@n,{k}],{k,1,n-1}]; (*sphere*)
AbsoluteTiming[chCx[bases]][[1]]
AbsoluteTiming[chCx0[bases]][[1]]
AbsoluteTiming[chCx1[bases]][[1]]
AbsoluteTiming[chCx2[bases]][[1]]


on a 2-core i3-560 CPU, I get:

0.774023

338.296

32.9942

17.921

Why is the nonparallelized code still the fastest? How can I efficiently parallelize this construction?

• chCx0 is incredibly slow because you're using a shared variable that gets updated very frequently. This makes it necessary for the parallel kernels to constantly communicate with each other, which has a lot of overhead. Shared variables should be used with caution. – Sjoerd Smit Nov 1 '18 at 21:22

There are several reasons for that, most notably the extensive use of AddTo (+=) with SparseArrays which makes parallelization through ParallelDo and friends impossible: Every AddTo will require a complete recomputation of the whole internal structure of the sparse arrays, so one process has to wait for the other in order to get write access.

Compared to the parallelism provided by linear algebra and vectorization, the Parallel-facilities are usually several orders of magnitude slower in constructing matrices.

Here is an implementation that shifts a lot of workload to vectorized routines (most notably using b.A to compute the facets of a simplex) and to the more efficient AssociationThread and Lookup. But most importantly, this constructs the combinatorics of the sparse arrays first and assembles them in one go.

getBoundaryMap[bases_] := Module[{dims, basesk, baseskk, bdr, m, b, A},
dims = Length /@ bases;
Table[
m = dims[[k]];
b = bases[[k]];
A = Normal[SparseArray[
Transpose[{
Flatten[Table[Delete[Range[k], i], {i, 1, k}]],
Range[1, k (k - 1)]
}] -> 1,
{k, k (k - 1)}]];
bdr = SparseArray[
Rule[
Transpose[{
Lookup[basesk, ArrayReshape[b.A, {m k, k - 1}]],
Flatten[Transpose[ConstantArray[Range[m], k]]]
}],
Flatten[ConstantArray[(-1)^Range[0, k - 1], m]]
],
{dims[[k - 1]], m}];
basesk = baseskk;
bdr,
{k, 2, Length[bases]}]
];


Here is a timing example for the case n = 16:

n = 16;
bases = Table[DeveloperToPackedArray@Subsets[Range@n, {k}], {k, 1, n - 1}];(*sphere*)
a = chCx[bases]; // AbsoluteTiming // First
b = getBoundaryMap[bases]; // AbsoluteTiming // First
a == b


37.8971

0.614126

True

Notice also the use of DeveloperToPackedArray; this is recommended because Subsets produces unpacked arrays and many routines will profit from packing them.

# Matrix trick

The main trick is to construct the matrix A for given k that such that b.A for an index vector $$b = (i_1,\dotsc,i_k)$$ returns the list

$$(\hat i_1,i_2,i_3, i_4, \dotsc,i_k, \quad i_1, \hat i_2, i_3, i_4,\dotsc,i_k, \quad i_1,i_2,\hat i_3, i_4,\dotsc,i_k, \quad \dotsc),$$

where $$\hat i_j$$ denotes an omitted index. Then it is only a matter of partitioning to obtain the list of facets:

$$( (\hat i_1,i_2,i_3, i_4, \dotsc,i_k), \quad (i_1, \hat i_2, i_3, i_4,\dotsc,i_k), \quad (i_1,i_2,\hat i_3, i_4,\dotsc,i_k), \quad \dotsc).$$

Here an example for $$k = 5$$:

k = 5;
b = Array[i, k];
A = Normal[
SparseArray[
Transpose[{Flatten[Table[Delete[Range[k], i], {i, 1, k}]],
Range[1, k (k - 1)]}] -> 1, {k, k (k - 1)}]];
ArrayReshape[b.A, {k, k - 1}] // MatrixForm


$$\left( \begin{array}{cccc} i(2) & i(3) & i(4) & i(5) \\ i(1) & i(3) & i(4) & i(5) \\ i(1) & i(2) & i(4) & i(5) \\ i(1) & i(2) & i(3) & i(5) \\ i(1) & i(2) & i(3) & i(4) \\ \end{array} \right)$$

The nice thing about this: If b is a list of index lists, b.A will compute all facets at once in vectorized way; in the end, it is again a matter of partitioning the subarrays with ArrayReshape:

ArrayReshape[b.A, {m k, k - 1}]

• Your method is very fast! On my laptop with i7-8850H CPU, case $n\!=\!23$, in which the largest matrix is bigger than $10^6\!\times\!10^6$, finishes in 102s, whilst my old method requires a few hours. But I don't understand the theory behind it, what is $A$? Would you be so kind and write something about this? Can this method be generalized to any chain complex, where columns (resp. rows) are indexed by basis elements $v$ (resp. $u$) via Association, and we are given a function $\partial(v)=\sum_u\alpha_{uv} u$ so that the matrix has $(u,v)$-th entry equal to $\alpha_{uv}$? – Leon Nov 2 '18 at 18:03
• @Leon See my edit. And yes, in principle, this should also work for other complexes. However, it will be important that bases already contains all possible simplices of the complex. – Henrik Schumacher Nov 2 '18 at 19:45
• Thank you very much! Give me a few days to absorb this... Also, your method seems to be non-parallelized, am I right? Why is chCx2 so slow? There are no SparseArray calculations, except once at the end. Are definitions not distributed properly? Should I have used SetSharedVariable instead? – Leon Nov 4 '18 at 0:39
• Hm. I guess one of the problems with chCx2  is the frequent call to Delete which involves a lot of fine grained memory allocation and copying. A second problem is that you manipulate lists of rules which cannot be done with packed array. Moreover, the calls to baseskk[s] are totally superfluous and the many calls to basesk[] are better done at once with Lookup. – Henrik Schumacher Nov 4 '18 at 7:47
• Associations are rather complicated data structure relying on external libraries. For all such data structures, it is a good strategy in Mathematica to reduce calls to them to a very least by calling them with as much data at a time as possible. Finally, I experience only few occasions where the parallelism of Parellel and friends was actually helpful for numerical computations. Usually, the parallelism provided by Compile and vectorized operations is much more efficient. That usually requires a bit different thinking and more coding effort but IMHO, it is totally worth it. – Henrik Schumacher Nov 4 '18 at 7:57