Parallelize the construction of sparse matrices

Question

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?

Answer 1

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;
   basesk = AssociationThread[bases[[1]], Range[dims[[1]]]];
   Table[
    m = dims[[k]];
    b = bases[[k]];
    baseskk = AssociationThread[b, Range[m]];
    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[Developer`ToPackedArray@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 Developer`ToPackedArray; 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}]