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

Question

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[123];
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.

Answer 1

Here is one idea using Pick and Total. Setup:

SeedRandom[123];
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.

Addendum

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}}

Answer 2

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]}.