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:
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 theGraphDistanceMatrix
part. Is there a way either the former or latter could be computed more efficiently?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.
a
and who ism
in your actual code (you could definitely make it clearer), but how does the naive approach usingPosition
andExtract
(e.g. something likeMean@Extract[m, Position[a, #, -1]] & /@ (DeleteDuplicates[ Flatten[a]] /. 0 -> Nothing)
) compare to yours in terms of timing? $\endgroup$distmatrix
in the code and $M$ ism
in the code. $\endgroup$