5
$\begingroup$

I have three lists $$A=\{\alpha_1...\alpha_N\}$$ $$B=\{\beta_1...\beta_N\}$$ $$M=\{M_1....M_N\}$$

In practice $\alpha_i$ and $\beta_i$ are complex numbers with modulus smaller than 1. The $M_i$s are integers in between -6 and 6.

I want to compute the following sum:

$$S=2\sum_i\sum_j\alpha_i\beta_j\delta_{M^i,M^j}$$

My best attempt is basically

ComputeSum[A_, B_,M_] := 
 Block[{kroneckerList},
  k= Outer[KroneckerDelta, M, M];
  ak= A*kroneckerList;
  Total[ak.B]
  ]

Which I believe does the job (but in any case what I want to compute is the sum above).
The problem is that $N$ can be quite big, up to $3\times 10^6$. and the k= Outer[KroneckerDelta, M, M] line gets computationally expensive.

Do you have faster solutions ?

$\endgroup$
2
  • $\begingroup$ Have you tried using simply Sum[..., {i, n}, {j, n}]? $\endgroup$ May 10, 2022 at 15:54
  • $\begingroup$ yes, that's a very bad way of doing it $\endgroup$
    – DarkBulle
    May 10, 2022 at 16:22

1 Answer 1

10
$\begingroup$
ComputeSum[A_, B_, M_] := A . Outer[KroneckerDelta, M, M] . B;

(*A helper function that for additive assembly of `SparseArray`s (_Mathematica's_ default is first in, last out.) *)
Options[MySparseArray] = {"Background" -> 0.};
MySparseArray[X_, r_, f_ : Total] := 
  If[(Head[X] === Rule) && (X[[1]] === {}),
   X[[2]],
   With[{spopt = SystemOptions["SparseArrayOptions"]},
    Internal`WithLocalSettings[
     SetSystemOptions[
      "SparseArrayOptions" -> {"TreatRepeatedEntries" -> f}],
     SparseArray[X, r, OptionValue["Background"]],
     SetSystemOptions[spopt]]
    ]
   ];

ComputeSum2[A_, B_, M_, k_] := Dot[
   MySparseArray[Partition[M + k + 1, 1] -> A, {2 k + 1}],
   MySparseArray[Partition[M + k + 1, 1] -> B, {2 k + 1}]
   ];

ComputeSum3[A_, B_, M_] := Dot[
   Values[GroupBy[Transpose[{M, A}], First -> Last, Total]],
   Values[GroupBy[Transpose[{M, B}], First -> Last, Total]]
   ];


n = 10000;
A = RandomReal[{-1, 1}, n];
B = RandomReal[{-1, 1}, n];
k = 6;
M = RandomInteger[{-k, k}, n];

result = ComputeSum[A, B, M]; // AbsoluteTiming // First
result2 = ComputeSum2[A, B, M, k]; // 
  AbsoluteTiming // First
result3 = ComputeSum3[A, B, M]; // AbsoluteTiming // First
Abs[result - result2]
Abs[result - result3]
Abs[result - result3]

16.8646

0.002846

0.006937

2.27374*10^-12

9.09495*10^-13

Edit

The idea of the two implementations is the same. We want to compute $$\begin{aligned} \sum_{i=1}^{n} \sum_{j=1}^n \alpha_i \, \delta_{M_i,M_j} \, \beta_j &= \sum_{i=1}^{n} \sum_{j=1}^n \sum_{k=-6}^6 \alpha_i \, \delta_{M_i,k} \, \delta_{k,M_j} \, \beta_j \\ &= \sum_{k=-6}^6 \left( \sum_{i=1}^{n}\alpha_i \, \delta_{M_i,k} \right) \, \left( \sum_{j=1}^n \delta_{k,M_j} \, \beta_j \right) \\ & = u^T v, \end{aligned}$$ where $$ u_k = \sum_{i=1}^n \alpha_i \, \delta_{M_i,k} \qquad v_k = \sum_{j=1}^n \beta_j \, \delta_{M_j,k}. $$ The naive summation costs $O(n^2)$; but each of u and v can be computed in $O((2\,k +1) \, n)$ time. So the new algorithm has complexity $$ O(2\,(2\,k +1) \,n + (2\,k +1)) = O(2\,(2\,k +1)\,n). $$ So if the range of k is much smaller than n, then we can save quite many flops this way.

Hence we may use

MySparseArray[Partition[M + k + 1, 1] -> A, {2 k + 1}]

(where we have to add shift the integers in M to be all greater than 0) or

Values[GroupBy[Transpose[{M, A}], First -> Last, Total]]

to assemble the vector u. Likewise we can do it for v. And in the end we just have Dot u and v together to get the result.

$\endgroup$
7
  • 1
    $\begingroup$ It's always you @Henrik Schumacher, thanks. This is some next level mathematica programming. Do you have a blog or an online class where you teach this kind of technique ? Because I know it works but I don't understand most your syntax $\endgroup$
    – DarkBulle
    May 10, 2022 at 10:50
  • 2
    $\begingroup$ No online material other than what I post here on StackExchange. But I added the mathematical explaination of why it works. I hope it helps. $\endgroup$ May 10, 2022 at 15:43
  • $\begingroup$ I now have to compute a sum that looks like this @Henrik Schumacher : $S=2\sum_i\sum_j\alpha_i\beta_j\delta_{M^i,M^j}\delta_{M'_i,M'_j}$, where $M'=\{M'_1....M'_N\}$. My code looks like this, where I use the idea of your answer without the first aprt which I don't understand very well. Can I do much better you think ? $\endgroup$
    – DarkBulle
    May 15, 2022 at 12:47
  • 1
    $\begingroup$ Well, $\delta_{M_i,M_j} \, \delta_{M'_i,M'_j} = \delta_{(M_i,M'_i),(M_j,M'_j)}$, so you can simply use ComputeSum3[A, B, Transpose[{M, M1}]]. The SparseArray thing would work, too, if appropriately adapted. I think it's Module[{numberOfPairs, knew, Mnew}, numberOfPairs = (2 k + 1)^2; knew = (numberOfPairs - 1)/2; Mnew = (2 k + 1) (M + k) + M1 + k - knew; ComputeSum2[A, B, Mnew, knew] ]. $\endgroup$ May 15, 2022 at 13:27
  • 1
    $\begingroup$ The idea behind the SparseArray trick here is to convert the pairs $(M_i,M'_i)$ into a single integer by treating $M_i+k+1$ and $M'_i + k + 1$ as digits of numbers with base $2 k +1$. $\endgroup$ May 15, 2022 at 13:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.