# Speed up construction of large matrix with if statements

I am constructing a matrix of the form $$T_{i_aj_ak_al_a,i_bj_bk_bl_b} = T_{i_ai_b}\delta_{j_aj_b}\delta_{k_ak_b}\delta_{l_al_b} + T_{j_aj_b}\delta_{i_ai_b}\delta_{k_ak_b}\delta_{l_al_b} + T_{k_ak_b}\delta_{i_ai_b}\delta_{j_aj_b}\delta_{l_al_b} + T_{l_al_b}\delta_{i_ai_b}\delta_{j_aj_b}\delta_{k_ak_b}$$

with

$$T_{i_ai_b} = \frac{\hbar^2 (-1)^{i_a-i_b}}{2m_w \Delta w^2} \left\{ \begin{array}{ll} \dfrac{\pi^2}{3}, & i_a = i_b \\ \dfrac{2}{(i_a-i_b)^2}, & i_a \neq i_b \end{array} \right\}$$

I did this in Mathematica with the following code:

Twxyz =
SparseArray[
Flatten[
Table[
If[ia == ib,
\[Pi]^2/3,
2/(ia - ib)^2]*(hbar^2*(-1)^(ia - ib))/(2 mw *dw^2)*
KroneckerDelta[ja, jb]*KroneckerDelta[ka, kb]*
KroneckerDelta[la, lb]
+
If[ja == jb,
\[Pi]^2/3,
2/(ja - jb)^2]*(hbar^2*(-1)^(ja - jb))/(2 mx*dx^2)*
KroneckerDelta[ia, ib]*KroneckerDelta[ka, kb]*
KroneckerDelta[la, lb]
+
If[ka == kb,
\[Pi]^2/3,
2/(ka - kb)^2]*(hbar^2*(-1)^(ka - kb))/(2 my*dy^2)*
KroneckerDelta[ia, ib]*KroneckerDelta[ja, jb]*
KroneckerDelta[la, lb]
+
If[la == lb,
\[Pi]^2/3,
2/(la - lb)^2]*(hbar^2*(-1)^(la - lb))/(2 mz*dz^2)*
KroneckerDelta[ia, ib]*KroneckerDelta[ja, jb]*
KroneckerDelta[ka, kb]
,
{lb, npointsgrid},
{la, npointsgrid},
{kb, npointsgrid},
{ka, npointsgrid},
{jb, npointsgrid},
{ja, npointsgrid},
{ib, npointsgrid},
{ia, npointsgrid}]
, {{1, 3, 5, 7}, {2, 4, 6, 8}}]
];


This approach works ok. The problem is that it takes several hours to build the matrix due to the many if statements it must process for each matrix element. Any suggestions on how to speed up the matrix construction?

EDIT:

The table is numeric.

Typical values for variables are

hbar = 6.62607004*^-34/(2*Pi)
npointsgrid = 13
mx*dx^2 = 3.76547*^-49


In the end, the matrix is constructed such that $$m_w\Delta w^2 = m_x\Delta x^2 = m_y \Delta y^2 = m_z \Delta z^2$$

In[307]:= mw*dw^2
mx*dx^2
my*dy^2
mz*dz^2

Out[307]= 3.76547*10^-49

Out[308]= 3.76547*10^-49

Out[309]= 3.76547*10^-49

Out[310]= 3.76547*10^-49


EDIT 2:

Correction of LaTeX code: should be $$\hbar^2$$ in definition of T, not $$\hbar$$ as previously written.

• You can use Piecewise. Oct 18, 2018 at 8:49

You're looking for KroneckerProduct:

hbar = 6.62607004*^-34/(2*Pi);
npointsgrid = 2;
mwdw2 = 3.76547*^-49;

npointsgrid = 13;

T = Table[((hbar^2 (-1)^(ia - ib)) Piecewise[{{π^2/3, ia == ib}}, 2/(ia - ib)^2])/(
2 mwdw2), {ia, npointsgrid}, {ib, npointsgrid}];

kro = SparseArray@Array[KroneckerDelta, {npointsgrid, npointsgrid}];

Twxyztest =
KroneckerProduct @@ RotateLeft[{T, kro, kro, kro}, #] & /@ Range@4 //
Total; // AbsoluteTiming
(* {0.109369, Null} *)


The validity can be easily checked by choosing a smaller npointsgrid. (For npointsgrid = 5, the original code takes about 12 seconds to finish calculating. )

• Thanks, that works really well. I knew there must be some product I could use because the matrix is rather sparse. Oct 19, 2018 at 12:22
• @de_thom_22 Glad it helps. If this answer solves your problem, you may accept it by clicking the checkmark sign. Oct 19, 2018 at 13:57