I have a list bdrs
of SparseArray
s, which represent the boundary matrices of a chain complex $R^{d_0}\overset{\partial_1}{\longleftarrow}R^{d_1}\longleftarrow\ldots\longleftarrow R^{d_{N-1}}\overset{\partial_N}{\longleftarrow}R^{d_N}$ over $R\!\in\!\{\mathbb{Z}, \mathbb{Q}, \mathbb{Z}_p\}$. I store this data as $d_0,\ldots,d_N$=dims
, $N$=dim
, $R$=p
$\!\in\!\{$"Z",0,p
$\}$. Typically, $\partial_k$ has dimensions $10^6\!\times\!10^6$ and density $10^{-5}$. See this or this or this link, to understand how SparseArray stores data.
For $k=1,\ldots,N$, I wish to (a) compute a sequence
Mm
=$M_1,\ldots,M_N$, in which $M_k$ is the set of all invertible entries in $\partial_k$, that have zeros left of them (in their row) and below them (in their column). Then for every $(i,j)$-entry in $M_k$, I wish to (b) delete it from $\partial_k$, and delete the $i$-column from $\partial_{k-1}$ and $j$-row from $\partial_{k+1}$.
First, if $R\!=\!\mathbb{Z}_p$ I mod out the entries, then delete all zero entries, and sort the rest of the entries:
If[p!="Z" && p!=0, bdrs=Mod[bdrs,p]];
bdrs = Table[SparseArray[Sort@ArrayRules@bdrs[[k]],dims[[k;;k+1]]], {k,dim}];
Attempt 1:
Note that some $d_k$ may be $0$. The construction is done by
Mm = DeleteCases[ Table[
i=bdrs[[k,All,j]]["NonzeroPositions"][[-1,1]]; w=bdrs[[k,i,j]];
If[bdrs[[k,i]]["NonzeroPositions"][[1,1]]==j && If[p=="Z", w^2==1, True],
bdrs[[k,i,j]]=0; {i,j}->w, ""],
{k,dim}, {j, If[dims[[k]]*dims[[k+1]]==0, {},
DeleteDuplicates@Flatten@bdrs[[k]]["ColumnIndices"]]}], "",2];
or
Mm = DeleteCases[ Table[ {i,j}=e; w=bdrs[[k,i,j]];
If[bdrs[[k,All,j]]["NonzeroPositions"][[-1]]=={i} && If[p=="Z", w^2==1, True],
bdrs[[k,i,j]]=0; {i,j}->w, ""],
{k,dim}, {e, If[dims[[k]]*dims[[k+1]]==0, {},
First/@ GatherBy[bdrs[[k]]["NonzeroPositions"], First]]}], "",2];
and then deletion of rows/columns:
Do[ If[1<k, bdrs[[k-1, All, #[[1,1]]&/@Mm[[k]] ]] = 0];
If[k<dim, bdrs[[k+1, #[[1,2]]&/@Mm[[k]], All ]] = 0];, {k,dim}];
Both constructions (a) are awfully slow on large matrices, but (b) is really fast.
Attempt 2:
Instead of a SparseArray
, I use an Association
of rows and of columns.
rows[bdr_] := With[
{l=GatherBy[ ArrayRules[bdr][[1;;-2]] /.({i_,j_}->w_Integer) :> {i,j,w}, First]},
Association@Table[r[[1,1]]->(Rest/@r), {r,l}]];
mm[r_,c_] := Module[{j,w}, Association@ DeleteCases[ Table[ {j,w}=r[i][[1]];
If[c[j][[-1]]=={i,w} && If[p=="Z",w^2==1,True], {i,j}->w,""],{i,Keys@r}],""]];
Mm = Table[r=rows@bdrs[[k]]; c=rows@Transpose@bdrs[[k]]; mm[r,c], {k,dim}];
Here (a) is much faster, but if I replace each bdrs[[k]]
with its association of rows and of columns, then achieving (b) becomes very slow (applying DeleteCases[#,MemberQ[#[[1,1]]&/@Mm[[k]],#]&]&
to all values in the Association
).
Comment:
Later, I will be computing a lot sums of rows, so I need to keep either bdrs
or rows & columns of these matrices as Association
s.
Examples for testing purposes:
Let us use a command, that builds up a chain complex from the faces of a simplicial complex:
chainComplexSC[bases_] := Module[{dim=Length@bases, dims=Length/@bases, basesk, baseskk, bdrs={}, entries, bdr, x},
basesk = AssociationThread[bases[[1]],Range@dims[[1]]];
Do[ baseskk=AssociationThread[bases[[k]],Range@dims[[k]]];
entries=Flatten[Table[{basesk[Delete[s,{{i}}]],baseskk[s]}->(-1)^(i+1),{s,Keys@baseskk},{i,k}],1];
AppendTo[bdrs, SparseArray[entries, dims[[k-1;;k]]]];
basesk=baseskk; Clear[baseskk];, {k,2,dim}]; bdrs];
sCxSimplices[facets_,k_] := ParallelCombine[ DeleteDuplicates@
Flatten[Table[Subsets[s,{k}],{s,#}],1]&, facets,Union,Method->"CoarsestGrained"];
Let us create the $m\!\times\!n$ chessboard complex:
m = 4;
n = 4;
facets = FindClique[GraphComplement@LineGraph@CompleteGraph[{m, n}], Infinity, All];
dim = Max[Length /@ facets];
bases = Table[sCxSimplices[facets, k + 1], {k, 0, dim - 1}];
dims = Length /@ bases;
bdrs = chainComplexSC[bases];
A good testing example is $m\!=\!8, n\!=\!9$, which has the f-vector dims
=$72, 2016, 28224, 211680, 846720, 1693440, 1451520, 362880$.