# Fast double sum involving Kronecker symbol

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 ?

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

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

(*A helper function that for additive assembly of SparseArrays (_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"]},
InternalWithLocalSettings[
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.

• 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 Commented May 10, 2022 at 10:50
• 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. Commented May 10, 2022 at 15:43
• 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 ? Commented May 15, 2022 at 12:47
• 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] ]. Commented May 15, 2022 at 13:27
• 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$. Commented May 15, 2022 at 13:31