2
$\begingroup$

I am trying something similar to the following code with NN as large as a few hundred (at least 100). Now it's very slow and most time is spent on calculating the defined summation function calling the eigenvalues and eigenvectors. Generating the eigensystem seems very fast.

But for such a simple form of function, I think it will still be very fast on C++ or so. Is there any way to speed up or overcome the bottleneck?

time0 = AbsoluteTime[];
NN = 40;
eigenR = Transpose[
   ParallelTable[
    Eigensystem[RandomReal[1, {2 NN, 2 NN}]], {k, NN}], {1, 3, 2}];
eigenL = Transpose[
  ParallelTable[
   Eigensystem[RandomReal[1, {2 NN, 2 NN}]], {k, NN}], {1, 3, 2}]; 
AbsoluteTime[] - time0

fun[nx1_, ny1_, nx2_, ny2_, e1_, e2_] := 
  Sum[(eigenR[[i, j, 1]])^-1 (eigenR[[i, j, 2, 2 nx1 + 1 + e1]] - 
      eigenR[[i, j, 2, 2 nx2 + 1 + e2]] Exp[
        I i (ny2 - ny1)]) Conjugate@(eigenL[[i, j, 2, 
        2 nx1 + 1 + e1]] - 
       eigenL[[i, j, 2, 2 nx2 + 1 + e2]] Exp[I i (ny2 - ny1)]), {i, 
    NN}, {j, 1, 2 NN}];
fundata = 
  ParallelTable[
   fun[nx, 0, nx, ny, 0, 1], {nx, 0, NN - 1}, {ny, 0, NN - 1}];
AbsoluteTime[] - time0
$\endgroup$
3
$\begingroup$

The following does the job approximately 250 times as fast on my machine. These steps were crucial:

  • Using Dot instead of Sum. This way, we may take advantage of vectorization and fuse-multiply-add instructions of the processor.

  • Keeping the arrays packed. That's crucial for profiting from vectorization. This was a bit of fiddling since Eigensystem returns unpacked arrays: The reason is that the input matrices are arrays of machie precision reals only (which can be packed) but the return values are mixed arrays of machine precision real and complex numbers. Therefore, I added 0. I to convert everything into machine precision comples numbers.

  • Transposing the arrays containing the eigenvectors once such that the read access can be done in the first slot. This way, the data read at a time lies in a contiguous array in memory which improves the streaming rate.

  • Conjugate the the L eigentvectors only once in the beginning.

  • Multiplying the inverses of the R eigenvalues with the R eigenvectors once in the beginning. This saves quite a lot multiplications and allows us to use a single Dot for the summation.

  • A slight speedup can be obtained by using Subtract instead of - because Mathematica interprets a - b as Plus[a, (-1.) b]], leading to a superflous multiplication operation.

But this had the most striking effect:

  • Inserting a 1. into Exp[I Range[NN] (ny2 - ny1)]. This way, the Range is constructed as packed array of machine precision reals. With integer array x, Exp[I x] stays unevaluated leading to a an unpacked array. This ruined actually all other efforts.

This is how the code looks now:

NN = 40;

Ldata = Developer`ToPackedArray /@ 
   Transpose[(Eigensystem /@ RandomReal[1, {NN, 2 NN, 2 NN}]) + 
     0. I];
LevecsConj = Conjugate[Transpose[Ldata[[2]], {2, 3, 1}]];

Rdata = Developer`ToPackedArray /@ 
   Transpose[(Eigensystem /@ RandomReal[1, {NN, 2 NN, 2 NN}]) + 
     0. I];
RevecsOverRevals = Transpose[
   Divide[
    Developer`ToPackedArray[Rdata[[2]]],
    Developer`ToPackedArray[Rdata[[1]]]
    ],
   {2, 3, 1}];

(* only for comparison*)
eigenL = Transpose[Ldata, {3, 1, 2}];
eigenR = Transpose[Rdata, {3, 1, 2}];

Your method:

fun[nx1_, ny1_, nx2_, ny2_, e1_, e2_] := Sum[
   Times[
    (eigenR[[i, j, 1]])^-1 , (eigenR[[i, j, 2, 2 nx1 + 1 + e1]] - 
      eigenR[[i, j, 2, 2 nx2 + 1 + e2]] Exp[I i (ny2 - ny1)]),
    Conjugate@(eigenL[[i, j, 2, 2 nx1 + 1 + e1]] - 
       eigenL[[i, j, 2, 2 nx2 + 1 + e2]] Exp[I i (ny2 - ny1)])
    ]
   , {i, NN}, {j, 1, 2 NN}];

    fundata = ParallelTable[ fun[nx, 0, nx, ny, 0, 1], {nx, 0, NN - 1}, {ny, 0, NN - 1}]; // 
      AbsoluteTiming // First

19.5267

Improved method

gun[nx1_, ny1_, nx2_, ny2_, e1_, e2_] := Dot[
   Flatten@Subtract[
     RevecsOverRevals[[2 nx1 + 1 + e1]],
     RevecsOverRevals[[2 nx2 + 1 + e2]] Exp[
       Range[1., NN] (I (ny2 - ny1))]
     ],
   Flatten@Subtract[
     LevecsConj[[2 nx1 + 1 + e1]],
     LevecsConj[[2 nx2 + 1 + e2]] Exp[Range[1., NN] (I (ny1 - ny2))]
     ]
   ];
gundata = ParallelTable[ gun[nx, 0, nx, ny, 0, 1], {nx, 0, NN - 1}, {ny, 0, NN - 1}]; // 
  AbsoluteTiming // First

0.079329

Check for accuracy:

Max[Abs[fundata - gundata]]

7.15277*10^-14

Notice tha both implementations have at least complexity $O(N^4)$ which is huge. So I suggest to think about what you actually want to do and to find out if there is an algorithm with complexity lower than $O(N^4)$.

$\endgroup$
  • $\begingroup$ Thank you so much for this excellent answer. When NN=100, the speedup is more than 50 times! Previously it takes more than an hour on a desktop, now I can do it on my laptop within a few minutes. But I'm not sure about the mechanism of this magic. There seem to be three things, PackedArray, vector elements moved to the first slot, and using Total & Times to rewrite the function. Do you have any idea which might be the most important one? $\endgroup$ – xiaohuamao Apr 16 at 5:06
  • $\begingroup$ @xiaohuamao I finally found the handbrake. See the newest version of my post. I also added detailed explanations. In total, keeping arrays packed is the most crucial part as only this will allow for vectorization. $\endgroup$ – Henrik Schumacher Apr 16 at 6:57
  • $\begingroup$ Thanks again for the nice answer and the update! $\endgroup$ – xiaohuamao Apr 16 at 12:34

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.