Efficient way of averaging over elements of one matrix based on another

Mathematical formulation of the problem with an example:

Suppose we have the matrix $$A$$ given by

$$A=\begin{pmatrix} 0 & 1 & 5\\ 1 & 0 & 5\\ 5 & 5 & 0 \end{pmatrix},$$

and the matrix $$M$$ generally of form:

$$M=\begin{pmatrix} 0 & m_{12} & m_{13}\\ m_{21} & 0 & m_{23}\\ m_{31} & m_{32} & 0 \end{pmatrix},$$

and the problem is to compute the arithmetic average of elements of $$M$$ where-ever the entries of $$A$$ are equal. That is, for each unique value of $$d$$ of $$A$$, finding the positions $$(i,j)$$ in $$A$$ where are $$a_{ij}=d,$$ and averaging over the entries of $$M$$ at the so found positions $$(i,j).$$ So in the above example we have $$2$$ unique values in $$A,$$ thus $$2$$ averages to compute over $$M,$$ namely:

• $$\langle m_1 \rangle= \frac{1}{2} (m_{12}+m_{21}),$$

• $$\langle m_5 \rangle= \frac{1}{4} (m_{13}+m_{23}+m_{31}+m_{32}).$$

Attempted approaches:

Method A

The matrices $$A$$ in my case are GraphDistanceMatrix[g] matrices, for a given undirected graph of $$n$$ nodes, with $$n$$ typically $$\approx 5 \times10^4,$$ and $$M$$ is a matrix of positive reals with zero diagonal and symmetric ($$A$$ is also always symmetric). So I figure, best bet to optimize the problem is to come up with a linear algebra formulation of it, e.g. if we can express the $$\langle m_i\rangle$$'s in terms of matrix products/operations of some sorts.

One approach I tried is:

• For each distinct value in $$A,$$ denoted by $$d$$ (e.g. $$d=1$$ in above example), we map all elements of $$A$$ to zero that are not equal to $$d.$$
• Then we count the number of non-zero elements left (needed to complete the arithmetic mean), so denoted by $$n_d$$.
• Then we take the Hadamard product between the transformed $$A$$ and $$M,$$ sum over all elements of resulting matrix and divide by $$n_d.$$

Here's an example:

SeedRandom;
n = 3000;
nedges = 12000;
g = RandomGraph[{n, nedges}];

distmatrix = GraphDistanceMatrix[g]; (*This is our matrix A*)
(*computing all distances. lengthy step 1*)

uniquedists =
Drop[Union@Flatten@distmatrix, 1];(*finding all distinct distances*)

d = uniquedists[[Length@uniquedists - 2]];

(*picking a distance as example, ultimately, we want to \
repeat what follows for all values in uniquedists*)

(*replacing all elements of distmatrix that are not equal to d by 0, \
and normalizing. Lengthy step 2*)

reduceddists = (distmatrix /. x_ /; x != d -> 0)/d; // AbsoluteTiming

{6.72683, Null}

(*counting how many nonzero left, which is the number of elements of \
m we'll be summing in the average, so nd is the normalization of our \
average*)
nd = Total@Flatten@reduceddists

2425810

(*Creating a random matrix m as example, which will be our M matrix in the problem statement.*)

m = RandomReal[{0.1, 1}, {n, n}];
m = UpperTriangularize[m, 1] + Transpose[UpperTriangularize[m, 1]];

(*made symmetric and its diagonal set to zero, to match our \
definition of M above.*)
v = ConstantArray[1, n];
(*now computing the average over all positions in m where element d \
was found in distmatrix*)

md = reduceddists*m; // AbsoluteTiming
(*Hadamard product to collapse all other values of m to zero*)

md = (v.md.v)/nd; // AbsoluteTiming
(*summing all remaining elements and normalizing*)

{0.611202, Null} {0.010056, Null}

md

0.550361

Question:

1. The bottleneck of my current approach above lies mostly in the replacements of values part, i.e.: (distmatrix /. x_ /; x != d -> 0)/d;, next to the GraphDistanceMatrix part. Is there a way either the former or latter could be computed more efficiently?

2. The above was my attempt to map the problem to a linear algebra one before solving it in Mathematica, but given the problem statement at the start, any other suggestions of approaches that can potentially be more efficient are most welcome. I feel like this whole thing can be done much more simply and I'm missing something obvious.

• I am not sure that I follow who is a and who is m in your actual code (you could definitely make it clearer), but how does the naive approach using Position and Extract (e.g. something like Mean@Extract[m, Position[a, #, -1]] & /@ (DeleteDuplicates[ Flatten[a]] /. 0 -> Nothing)) compare to yours in terms of timing? Nov 20 '19 at 18:22
• @MarcoB you're right, I should be clearer about those definitions as I used different variable names than those in the statement. $A$ is distmatrix in the code and $M$ is m in the code.
– user52181
Nov 20 '19 at 18:23
• Thank you for the clarification! Nov 20 '19 at 18:26

Here is one idea using Pick and Total. Setup:

SeedRandom;
n = 3000;
nedges = 12000;
g = RandomGraph[{n,nedges}];

distmatrix=GraphDistanceMatrix[g];

m = RandomReal[{0.1,1}, {n, n}];
m = UpperTriangularize[m, 1]+Transpose[UpperTriangularize[m, 1]];
ones = ConstantArray[1, {n, n}];

distances = Union @ Flatten @ distmatrix

{0, 1, 2, 3, 4, 5, 6, 7}

Now, using Pick and Total to find the averages for 6:

Total[Pick[m, distmatrix, 6], Infinity]/Total[Pick[ones, distmatrix, 6], Infinity] //AbsoluteTiming

{0.328108, 0.550647}

A table of the values for the positive distances:

Table[Total[Pick[m, distmatrix, i], Infinity]/Total[Pick[ones, distmatrix, i], Infinity], {i, Rest @ distances}] //AbsoluteTiming

{1.93576, {0.548076, 0.548365, 0.549827, 0.549896, 0.550361, 0.550647, 0.527063}}

Another possibility is to create a pair matrix and use GroupBy. To be most efficient, though, the distances should be reals. So:

pairs = Flatten[{N @ distmatrix, m}, {{2,3}, {1}}]; //AbsoluteTiming

{0.305517, Null}

creates a list of {d, v} pairs. Then using GroupBy:

r2 = KeySort @ KeyMap[Round] @ GroupBy[pairs, First -> Last, Mean]; //AbsoluteTiming
Rest @ Values @ r2

{1.58506, Null}

{0.548076, 0.548365, 0.549827, 0.549896, 0.550361, 0.550647, 0.527063}

which is the same as before.

Note that the GroupBy approach can be sped up a bit further by making use of the "GroupByList" ResourceFunction:

ResourceFunction["GroupByList"][Flatten @ m, Flatten @ distmatrix, Mean] //AbsoluteTiming

{0.148771, <|0 -> 0., 4 -> 0.549896, 5 -> 0.550361, 3 -> 0.549827, 2 -> 0.548365, 6 -> 0.550647, 1 -> 0.548076, 7 -> 0.527063|>}

An even faster possibility is to work with vectors instead (modified from earlier answer by removing the unnecessary Unitize)

fd = Flatten @ distmatrix;
fm = Flatten @ m;

Table[
With[{v = Clip[fd, {i, i}, {0, 0}]},
v . fm / Total[v]
],
{i, Rest @ distances}
] //AbsoluteTiming

{0.445384, {0.548076, 0.548365, 0.549827, 0.549896, 0.550361, 0.550647, 0.527063}}

• Thanks a lot! It's incredible how effortless and simple you make it look Carl, these suggestions are simply brilliant! I have also tested these approaches on larger graphs, and your 2nd method using GroupBy seems to perform the fastest, about $3$ times faster (in absolute time) than the vectorial approach.
– user52181
Nov 21 '19 at 8:54
• @user929304 Try out the "GroupByList" resource function, it should be faster than the GroupBy approach. Nov 21 '19 at 17:18

One way to do this is to set up a matrix that is 1 where everything is the same and zero where different, then Total all the elements in the product and divide by how many there are to get the mean. For your example A and M:

a = {{0, 1, 5}, {1, 0, 5}, {5, 5, 0}};
M = Array[m, {3, 3}];
sel = Union[Flatten[a]];
Table[t = Total[(1 - Abs@Sign[a - sel[[i]]]) M, 2]; t/Length[t], {i, 1, Length[sel]}]

{1/3 (m[1, 1] + m[2, 2] + m[3, 3]), 1/2 (m[1, 2] + m[2, 1]),
1/4 (m[1, 3] + m[2, 3] + m[3, 1] + m[3, 2])}

This gives a list with the three outputs corresponding to the means of all the 0s, then the 1s, then the 5s. If you don't want the first term, you can change the iterator to {i, 2, Length[sel]}.

• Thanks Bill, this is somewhat close to what I tried first. I've tested your approach on the random graph used in OP, it seems to encounter division by zero.
– user52181
Nov 21 '19 at 8:59