Revisions to Fast double sum involving Kronecker symbol

deleted 65 characters in body

Source Link

edited May 15, 2022 at 13:14

109.4k
7
186
322

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, \[Alpha]A}], First -> Last, Total]],
   Values[GroupBy[Transpose[{M, \[Beta]B}], First -> Last, Total]]
   ];


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

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

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

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, \[Alpha]}], First -> Last, Total]],
   Values[GroupBy[Transpose[{M, \[Beta]}], First -> Last, Total]]
   ];


n = 10000;
\[Alpha] = RandomReal[{-1, 1}, n];
\[Beta] = RandomReal[{-1, 1}, n];
k = 6;
M = RandomInteger[{-k, k}, n];

result = ComputeSum[\[Alpha], \[Beta], M]; // AbsoluteTiming // First
result2 = ComputeSum2[\[Alpha], \[Beta], M, k]; // 
  AbsoluteTiming // First
result3 = ComputeSum3[\[Alpha], \[Beta], M]; // AbsoluteTiming // First
Abs[result - result2]
Abs[result - result3]
Abs[result - result3]

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

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]

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

edited body

Source Link

edited May 10, 2022 at 15:41

Henrik Schumacher

109.4k
7
186
322

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)$O(n^2)$; but each of u and v can be computed in O((2\,k +1) \, n)$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.

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.

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.

\; added 293 characters in body; added 8 characters in body

Source Link

edited May 10, 2022 at 11:27

Henrik Schumacher

109.4k
7
186
322

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, \[Alpha]}], First -> Last, Total]],
   Values[GroupBy[Transpose[{M, \[Beta]}], First -> Last, Total]]
   ];


n = 10000;
\[Alpha] = RandomReal[{-1, 1}, n];
\[Beta] = RandomReal[{-1, 1}, n];
k = 6;
M = RandomInteger[{-k, k}, n];

result = ComputeSum[\[Alpha], \[Beta], M]; // AbsoluteTiming // First
result2 = ComputeSum2[\[Alpha], \[Beta], M, k]; // 
  AbsoluteTiming // First
result3 = ComputeSum3[\[Alpha], \[Beta], M]; // AbsoluteTiming // First
Abs[result - result2]
Abs[result - result3]
Abs[result - result3]

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, \[Alpha]}], 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.

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

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, \[Alpha]}], First -> Last, Total]],
   Values[GroupBy[Transpose[{M, \[Beta]}], First -> Last, Total]]
   ];


n = 10000;
\[Alpha] = RandomReal[{-1, 1}, n];
\[Beta] = RandomReal[{-1, 1}, n];
k = 6;
M = RandomInteger[{-k, k}, n];

result = ComputeSum[\[Alpha], \[Beta], M]; // AbsoluteTiming // First
result2 = ComputeSum2[\[Alpha], \[Beta], M, k]; // 
  AbsoluteTiming // First
result3 = ComputeSum3[\[Alpha], \[Beta], M]; // AbsoluteTiming // First
Abs[result - result2]
Abs[result - result3]
Abs[result - result3]

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, \[Alpha]}], First -> Last, Total]],
   Values[GroupBy[Transpose[{M, \[Beta]}], First -> Last, Total]]
   ];


n = 10000;
\[Alpha] = RandomReal[{-1, 1}, n];
\[Beta] = RandomReal[{-1, 1}, n];
k = 6;
M = RandomInteger[{-k, k}, n];

result = ComputeSum[\[Alpha], \[Beta], M]; // AbsoluteTiming // First
result2 = ComputeSum2[\[Alpha], \[Beta], M, k]; // 
  AbsoluteTiming // First
result3 = ComputeSum3[\[Alpha], \[Beta], M]; // AbsoluteTiming // First
Abs[result - result2]
Abs[result - result3]
Abs[result - result3]

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, \[Alpha]}], 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.

Source Link

created May 9, 2022 at 23:00

Henrik Schumacher

109.4k
7
186
322

Loading

Stack Exchange Network

Return to Answer