# Counting non-zero multiplications in large sparse matrix multiplications

Suppose I have several large, fixed, sparse matrices of real-valued entries. To make things simple, assume each is $$n \times n$$. I will multiply them, e.g., matA.matB.matC. Of course I can reduce the computational complexity significantly because of the sparse nature of the matrices, avoiding all multiplications by $$0$$.

My question is how to algorithmically count the number of such required multiplications for a product of (say) three given sparse matrices.

To illustrate the computational savings I can work with known matrices that are not SparseMatrix, and a typical value is $$n = 32$$.

My first approach was to write out the multiplications using summations over indexes, and such, but it started getting involved and complicated. I can use Position[myA, Except[0],2] to get the locations of non-zero entries in the matrix. But then what?

I thought there might be a more elegant solution.

• With sparse matrices, Position[myA, Except[0],2] doesn't work, but you can find those positions with myA["NonzeroPositions"]. Jan 18 at 23:58
• maybe Total[Unitize[a].Unitize[b], 2]?
– kglr
Jan 19 at 13:23
• ... and for three matrices Total[Fold[Unitize[#].#2&,Unitize/@{a,b,c}], 2]?
– kglr
Jan 19 at 21:51
• Oh yes... that works! (Another elegant solution from the incomparable @kglr.). Care to post this so we can advertise your solution and close further discussion? Jan 19 at 21:55
• Oh I see. I thought you wanted to count the flops without performing the flops (or the according number of integer operations). Jan 20 at 1:48

$$\{i,j\}^{th}$$ entry in a.b is $$a_{\text{i}1} b_{1\text{j}}+a_{\text{i}2} b_{2\text{j}}+a_{\text{i}3} b_{3\text{j}}+\cdots + a_{\text{i}n} b_{n\text{j}}$$.

So, if we Unitize the input matrices, for ua = Unitize[a] and ub = Unitize[b], $$\{i,j\}^{th}$$ entry in ua.ub gives the number of terms where both multiplicands are non-zero. Therefore, to get the total number of multiplications we add the entries in ua.ub:

Total[Unitize[a].Unitize[b], 2]


For three matrices,a, b and c, we can use this method iteratively:

Total[#, 2] & /@ Rest[FoldList[Unitize[#] . #2 &, Unitize /@ {a, b, c}]]


to get the number of multiplications in ab = a.b and ab.c.

• (+1) Nicely done, @kglr! Jan 19 at 23:50
• Thank you @E.Chan-López.
– kglr
Jan 19 at 23:59
• Yep. ($\checkmark$). Jan 20 at 4:41

Say your matrices are a and b. You have a nonzero multiplication for each pair of entries {a[[i,k]],b[[k,j]]} in which both a[[i,k]] and b[[k,j]] are nonzero. So for each k we merely have to count the number of nonzero entries in column k of a and the number of entries in row k of b and then multiply these counts. Afterwards, we some over all k.

Getting the number of nonzeroes for each row of b is easily optained by

Differences[b["RowPointers"]]


To get the number of nonzeros per column of a we run

Differences[Transpose[a]["RowPointers"]]


In order to multiply and sum these counts, we just need to do this:

Dienter code herefferences[Transpose[a["RowPointers"]].Differences[b["RowPointers"]]


If you don't want to transpose the matrix a, then you can use the following algorithm instead: You go through each row i of a; for each nonzero entry at position {i,k} you add the number of nonzeros of row k of b to a counter. Here is an implementation of it:

cCountNonzeroMultiplication =
Compile[{{arp, _Integer, 1}, {aci, _Integer, 2}, {brp, _Integer, 1}},
Module[{counter, m, k},
m = Length[arp] - 1;
counter = 0;
Do[
Do[
k = CompileGetElement[aci, l, 1];
counter +=
CompileGetElement[brp, k + 1] - CompileGetElement[brp, k]
, {l, CompileGetElement[arp, i] + 1,
CompileGetElement[arp, i + 1]}]
, {i, 1, m}];
counter
],
CompilationTarget -> "C",
RuntimeOptions -> "Speed"
];


This gives you the flop count for the product of two matrices. For this we need either

• the rowpointers of one transposed matrix and the rowpointers of the other
• the rowpointers of both matrices.

I don't know, yet, how to do that for three matrices a, b, c without computing the sparsity pattern of a.b and b.c first. In principle, I know how to do compute the sparsity pattern. But that requires quite precisely as many integer operations as a.b requires flops. =/

• Hi, @Henrik Schumacher! I suppose it's possible to build an iterative process from your answer to do the counting when dealing with the product of many matrices. Jan 19 at 6:19
• Hm. I don't know. See the remark I added to my post. Jan 19 at 12:10
• Thanks ($+1$), and an excellent step in the right direction. Jan 19 at 17:20

To get the counts without performing multiplications, we can also use the SparseArray property "NonzeroPositions" as follows:

$$\{i,j\}^{th}$$ entry in a.b:

$$a_{\text{i}1} b_{1\text{j}}+a_{\text{i}2} b_{2\text{j}}+a_{\text{i}3} b_{3\text{j}}+\cdots + a_{\text{i}n} b_{n\text{j}}$$.

The pair $$a_{ik}$$ and $$b_{kj}$$ are both non-zero iff the position {i,k} belongs to a["NonzeroPositions"] and {k,j} belongs to b["NonzeroPositions"].

So we can take tuples of non-zero positions of a and b and count the pattern {{_,k_},{k_,_}} to get the desired number of multiplications with non-zero multiplicands in a.b:

{nzpa, nzpb, nzpc} = #["NonzeroPositions"] & /@ {a, b, c}

Count[{{_, k_}, {k_, _}}] @ Tuples[{nzpa, nzpb}]


Non-zero positions in ab = a.b can be obtained as

nzpadotb = DeleteDuplicates @ Cases[{{i_, k_}, {k_, j_}} :> {i, j}] @
Tuples[{nzpa, nzpb}]


So, to get the number of multiplications with non-zero multiplicands in ab . c is

Count[{{_, k_}, {k_, _}}] @ Tuples[{nzpadotb, nzpc}]


Note: In recent versions, you can also use the property name "ExplicitPositions" in place of "NonzeroPositions".

• Sorry, but I have to say this: Why do you post this totally inefficient, memory-wasting piece of code after you posted the other really good answer that simply uses efficient sparse matrix routines? Already for medium sized sparse arrays with about 10000 nonzero entries, the result of Tuples[{nzpa, nzpb}]` just does not fit into memory! Jan 20 at 13:06
• @HenrikSchumacher, excellent point as always -- thank you. I simply missed "... large ..." in the very first line of David's question and assumed n =32 was the size of input matrices. Thought the challenge in your comment (get the desired result using only the properties of input matrices without taking matrix dot products) was too interesting not to take up. ... and, if mmase had tracked the number of "totally inefficient, memory-wasting piece of code" posted by users, I probably would have ranked really high -- if not the highest:)
– kglr
Jan 20 at 23:15