Table Creation Slow Performance

Question

The following is a question in regards to the matrix $A$ (originally the zero $4n \times 4n$ matrix) below:

 In[288]:= For[j = 1, j <= 2 n, j++,
   For[k = 1, k <= 2 n, k++,
     A[[2 j - 1]][[2 k - 1]] = -2 I*H[[j]][[k]] - M[[k]][[j]] + 
       M[[j]][[k]];
     A[[2 j - 1]][[2 k]] = 2*I*M[[j]][[k]];
     A[[2 j]][[2 k - 1]] = -2*I*M[[k]][[j]];
     A[[2 j]][[2 k]] = -2 I*H[[j]][[k]] + M[[k]][[j]] - M[[j]][[k]];
     ];
   ]; // AbsoluteTiming

Out[288]= {25.767, Null}

$H$ and $M$ are both matrices of size $2n \times 2n$ and $n = 640$.

I have tried using Table (code below) but it seems to run much slower than the above code.

In[291]:= A = Table[
    p = Ceiling[j/2]; q = Ceiling[k/2];

    If[OddQ[j] && OddQ[k], -2 I*H[[p]][[q]] - M[[q]][[p]] + 
       M[[p]][[q]], 0]
     If[OddQ[j] && EvenQ[k], 2*I*M[[p]][[q]], 0]
     If[EvenQ[j] && OddQ[k], -2*I*M[[q]][[p]], 0]
     If[EvenQ[j] && EvenQ[k], -2 I*H[[p]][[q]] + M[[q]][[p]] - 
       M[[p]][[q]], 0]

    , {j, 1, 4 n}, {k, 1, 4 n}]; // AbsoluteTiming

Out[291]= {52.4175, Null}

I was wondering what possible speed ups there are for the generation of $A$?

Answer 1

n = 640;
A1 = A2 = A3 = ConstantArray[0. + 0. I, {4 n, 4 n}];
{H, M} = RandomComplex[{-1 - I, 1 + I}, {2, 2 n, 2 n}];

The original code and its timing from my machine:

For[j = 1, j <= 2 n, j++,
    For[k = 1, k <= 2 n, k++,
      A1[[2 j - 1]][[2 k - 1]] = -2 I*H[[j]][[k]] - M[[k]][[j]] + 
        M[[j]][[k]];
      A1[[2 j - 1]][[2 k]] = 2*I*M[[j]][[k]];
      A1[[2 j]][[2 k - 1]] = -2*I*M[[k]][[j]];
      A1[[2 j]][[2 k]] = -2 I*H[[j]][[k]] + M[[k]][[j]] - 
        M[[j]][[k]];
      ];
    ]; // AbsoluteTiming // First

39.3159

The same with a Do loop and with simplified Part-syntax (main speedup is indeed achieved by the change of the Part-syntax!):

Do[
    A2[[2 j - 1, 2 k - 1]] = -2 I H[[j, k]] - M[[k, j]] + M[[j, k]];
    A2[[2 j - 1, 2 k]] = 2 I M[[j, k]];
    A2[[2 j, 2 k - 1]] = -2 I M[[k, j]];
    A2[[2 j, 2 k]] = -2 I H[[j, k]] + M[[k, j]] - M[[j, k]];
    , {j, 1, 2 n}, {k, 1, 2 n}]; // AbsoluteTiming // First

18.3622

And a vectorized approach:

First@RepeatedTiming[
  A3[[1 ;; ;; 2, 1 ;; ;; 2]] = -2 I H - Transpose[M] + M;
  A3[[1 ;; ;; 2, 2 ;; ;; 2]] = 2 I M;
  A3[[2 ;; ;; 2, 1 ;; ;; 2]] = -2 I Transpose[M];
  A3[[2 ;; ;; 2, 2 ;; ;; 2]] = -2 I H + Transpose[M] - M;
  ]
A1 == A2 == A3

0.150

For the meaning of ;;, lookup Span in the documentation.

Also notice that it is important to initialize A3 as matrix of complex zeroes 0. + 0. I in order to keep A3 a packed array.