3
$\begingroup$

I created two loops using Parallel table and a table. The outer loop causes a (nstates by nstates) dimension matrix and the inner loop calculates the matrix elements based on their position in the outer matrix and also the initial predefined array (avec). Keeping in mind the avec I wrote here is just an example and always not so simple, and I do not usually end up with a final diagonal matrix, I want to know how I can make use of SparsArray and probably Associate to speed it up? I have tried to use SparseArray, but since I need to calculate the elements initially based on their position and I am not that familiar with the tricks to do that, I would appreciate it if someone can help.

here is the example code:

initial conditions:

ℓ0 = 8;
γ = 
Join[Table[{m, 1}, {m, -ℓ0, ℓ0}], 
Table[{m, -1}, {m, -ℓ0, ℓ0}]]
nstates = Length[γ]
ne = 2 ℓ0 + 2

simple avec:

 avec = Table[0, {ie, 1, ne}, {i, 1, nstates}];
 Do[avec[[1, 3]] = 1;
 avec[[i + 1, 2 + ne]] = 1, {i, 1, (ne - 1)}]

This is the function I am using in the loop, and probably its at its best speed, so I do not think I need to change anything about this part:

 ParallelEvaluate[Off[ClebschGordan::phy];
 ClearAll[j3s];
 j3s[a_, b_, c_] := j3s[a, b, c] = ThreeJSymbol[a, b, c];
 ClearAll[dfxn]; 
 dfxn[ℓ_, m1_, m2_, p1_, p2_] := 
 N@If[m1 + p1 == m2 + p2, 
 Sum[(2 ℓ + 1)^2 (2 ℓtemp + 1)/(4 π )
    Sum[If[m1 + p1 == mval && m2 + p2 == mval, 
     j3s[{ℓ, m1}, {ℓ, 
        p1}, {ℓtemp, -mval}] j3s[{ℓ, 
        m2}, {ℓ, 
        p2}, {ℓtemp, -mval}] j3s[{ℓ, 
        0}, {ℓ, 0}, {ℓtemp, 0}]^2, 
     0], {mval, -ℓtemp, ℓtemp}], \
 {ℓtemp, 0, 2 ℓ}], 0];];

and here is the loop which I want to modify the speed using SparseArray if possible:

   vex =(*(2 ℓ0 +1)^2*) ParallelTable[

  mpf = γ[[f, 1]];
  mk = γ[[k, 1]];
  μpf = γ[[f, 2]];
  μk = γ[[k, 2]];
  Chop[Total[
  Table[(* Here we loop over the HF states 
  *)Off[ClebschGordan::phy];
  pi = γ[[i, 1]]; 
  pj = γ[[j, 1]];
  μpi = γ[[i, 2]];
  μpj = γ[[j, 2]];
  If[μpi == μk && μpj == μpf, 
  N[Conjugate[avec[[ie, i]]] *avec[[ie, j]]* 
    dfxn[ ℓ0, pi, mk, mpf, pj]], 0]
 , {ie, 1, ne}, {i, 1, nstates}, {j, 1, nstates}]
 , Infinity]],
 {f, 1, nstates}, {k, 1, nstates}]
------------------------------------------------------------------

This was the simple version of my program, Thanks to you @Henrik Schumacher I got to learn about packed arrays. But I am still having problem with making it work for different avecs in which I may have more than one nonzero element in each row. for example if I define my avec as

 avec = SparseArray@ConstantArray[0., {ne, nstates}];

 stateList = Flatten[Table[stateA1 = {im - 1, -1}; stateA2 = {-im, 1};
 stateB1 = {-im, -1}; stateB2 = {im, 1};
 iA1 = Part[Position[\[Gamma], stateA1], 1, 1];
 iA2 = Part[Position[\[Gamma], stateA2], 1, 1];
 iB1 = Part[Position[\[Gamma], stateB1], 1, 1];
 iB2 = Part[Position[\[Gamma], stateB2], 1, 1];
 {{iA1, iA2}, {iB1, iB2}}, {im, 1, \[ScriptL]0}], 1];
 Do[avec[[ie, stateList[[ie, 1]]]] = 
 Sin[(ie \[Pi])/(2 (2 \[ScriptL]0 + 1))];
 avec[[ie, stateList[[ie, 2]]]] = 
 Cos[(ie \[Pi])/(2 (2 \[ScriptL]0 + 1))];
, {ie, 1, 2 \[ScriptL]0}]
 avec[[2 \[ScriptL]0 + 2, Part[Position[\[Gamma], {0, 1}], 1, 1]]] = 1;
 avec[[2 \[ScriptL]0 + 1, 
 Part[Position[\[Gamma], {\[ScriptL]0, -1}], 1, 1]]] = 1;

Then I don't know how to change this part of your code

 aa = ConjugateTranspose[avec].avec;
 {ilist, jlist} = Transpose[aa["NonzeroPositions"]];

Or if anyway I can use the same method still. Thank you again for your time and help.

$\endgroup$
  • 1
    $\begingroup$ Yes, it is certainly possible to speed this up considerably. Basically 99.9 percent of the numbers in the Table you Total over are zeroes. There should be a simple logic to figure out whether the entries have to be computed in the first place. To point out what I mean: constructions like Sum[f[i] g[j] KroneckerDelta[i, j], {i, 1, n}, {j, 1, n}] are popular with physicists as they allow simple paper calculation, but super inefficient on a compute because they require $\Theta(n^2)$. Instead, Sum[f[i] g[i], {i, 1, n}] leads to the same result in $\Theta(n)$. $\endgroup$ – Henrik Schumacher Mar 17 '20 at 6:59
  • 1
    $\begingroup$ Another issue (related to that) is that the final matrix seems to be a diagonal matrix. So there is also no need to compute the off-diagonal entries. $\endgroup$ – Henrik Schumacher Mar 17 '20 at 6:59
5
$\begingroup$

Okay, here is my take on it, but I give no guarantees for correctness. The basic idea here is to replace summations by logic, to do as much logic as possible offline and to do remaining by filtering lists of summation indices.

ℓ0 = 8;
γ = Developer`ToPackedArray@ Join[Table[{m, 1}, {m, -ℓ0, ℓ0}], Table[{m, -1}, {m, -ℓ0, ℓ0}]];
{γ1, γ2} = Transpose[γ];
nstates = Length[γ];
ne = 2 ℓ0 + 2;

avec = SparseArray@ConstantArray[0., {ne, nstates}];
Do[avec[[1, 3]] = 1.; avec[[i + 1, 2 + ne]] = 1., {i, 1, (ne - 1)}];

aa = ConjugateTranspose[avec].avec;
{ilist, jlist} = Transpose[aa["NonzeroPositions"]];
vals = aa["NonzeroValues"];
γ1ilist = γ1[[ilist]];
γ1jlist = γ1[[jlist]];
γ2ilist = γ2[[ilist]];
γ2jlist = γ2[[jlist]];

ParallelEvaluate[Off[ClebschGordan::phy];
  ClearAll[j3s];
  j3s[a_, b_, c_] := j3s[a, b, c] = N@ThreeJSymbol[a, b, c];
  ClearAll[dfxn];
  dfxn[ℓ_, m1_, m2_, p1_, p2_] := 
   If[m1 + p1 == m2 + p2, 
    Sum[(2 ℓ + 1)^2 (2 ℓtemp + 1)/(4 π) Sum[
       If[m1 + p1 == mval && m2 + p2 == mval, 
        j3s[{ℓ, m1}, {ℓ, p1}, {ℓtemp, -mval}] j3s[{ℓ, m2}, {ℓ, p2}, {ℓtemp, -mval}] j3s[{ℓ, 0}, {ℓ, 0}, {ℓtemp, 0}]^2, 
        0.], {mval, -ℓtemp, ℓtemp}], 
{ℓtemp, 0, 2 ℓ}], 0.];];

And now

vex2 = SparseArray@ParallelTable[
    γ1f = γ1[[f]];
    γ1k = γ1[[k]];
    γ2f = γ2[[f]];
    γ2k = γ2[[k]];
    idx = Intersection[
      Random`Private`PositionsOf[γ2ilist, γ2k],
      Random`Private`PositionsOf[γ2jlist, γ2f]
      ];
    vals[[idx]].MapThread[
      dfxn[ℓ0, #1, γ1k, γ1f, #2] &,
      {γ1ilist[[idx]], γ1jlist[[idx]]}
      ],
    {f, 1, nstates}, {k, 1, nstates}];

This gets the job done in about 0.15 s (OP's code required about 100 s) on my machine. For the data provided by OP, vex and vex2 coincide. I am pretty sure that dfxn can be optimized further; it still contains summation of many zeroes. But I leave it as it is for the moment.

$\endgroup$
  • $\begingroup$ Thank you so much Henrik, I didn't know about packed-arrays. $\endgroup$ – Delaram Nematollahi Mar 19 '20 at 1:53
  • 1
    $\begingroup$ This is indeed odd. The problem seems to be that the values of avec are not all floating point number anymore; for some reason I do not know, this drives_Mathematica_ to the assumption that ConjugateTranspose[avec].avec is a dense matrix. This can easily circumvented with aa = ConjugateTranspose[N@avec].N@avec; or by enforcing that the values that you write into avec are all floating point numbers. $\endgroup$ – Henrik Schumacher Mar 23 '20 at 6:54
  • 1
    $\begingroup$ Yeah. That won't work. But if your condition is μpi == μpj && μk== μpf, then you do not have to cycle over j and f at all: Just put μpj = μpi and μpf = μk. This spares you two loops. $\endgroup$ – Henrik Schumacher Mar 24 '20 at 6:20
  • 2
    $\begingroup$ The overall point to efficiency is this: First find out which combinations of indices lead to a nonzero result -- without(!) touching all possible combinations. Then do the computations only for those indices. There are different strategies to do that, depending on the type of conditions. And of course, it is impossible for me to anticipate all of them. I am afraid, you will have to learn yourself how to figure that out in the long run. $\endgroup$ – Henrik Schumacher Mar 24 '20 at 6:20
  • 2
    $\begingroup$ Well, I think one has to learn them by doing. I don't know of any particular reference (which is not to say there weren't any). I suggest to do the following: Whenever you see a check for equality or nonzeroeness in your code, ask yourself whether you could have known the outcome offline. If you do, that gives you a hint how to avoid superfluous loops. And be greedy very greedy with nested loops; they let your runtime explode. $\endgroup$ – Henrik Schumacher Mar 24 '20 at 17:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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