# How to speed up my code using Sparse-array possibly?

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.

• 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)$. Mar 17, 2020 at 6:59
• 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. Mar 17, 2020 at 6:59

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;
γ = DeveloperToPackedArray@ 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[
RandomPrivatePositionsOf[γ2ilist, γ2k],
RandomPrivatePositionsOf[γ2jlist, γ2f]
];

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.
• 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. Mar 23, 2020 at 6:54
• 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. Mar 24, 2020 at 6:20