I have two lists. The first, coefficientlist, is dense and has dimension {176,176} elements. The second, cubelist, is a list of three-dimensional lists, i.e. cubelist has dimension {176,65,65,65}. I need to perform the following summation:

      sumlist += coefficientlist[[n,m]] (cubelist[[n]]cubelist[[m]]) ];

Essentially I want to perform a sum like $S = \sum\limits_{i ,j} T_{ij} \phi_i \phi_j $ where the $T$ matrix is called coefficientlist in the code above and the $\phi_i$ are represented as grids over real space.

I have to do this several times, so I want to find the most efficient way to do it.

Thanks in advance.

  • $\begingroup$ Your code simply doesn't work for me... Table[0,{65},{65},{65}] generates a large matrix of... zeros! So the sum must be zero too... $\endgroup$
    – Rod
    Jun 4 '13 at 18:19
  • 2
    $\begingroup$ @Rod Lm The sum is zero even if the matrix is really, really large? ;} $\endgroup$
    – DavidC
    Jun 4 '13 at 18:26
  • $\begingroup$ @RodLm I think sumlist is a container to record the sum, sumlist=Table[0,{65},{65},{65}] is merely an initialization. $\endgroup$
    – Silvia
    Jun 4 '13 at 18:27
  • 1
    $\begingroup$ Another question: if cubelist is three-dimensional, how can it have Dimension {176,65,65,65}? $\endgroup$
    – Rod
    Jun 4 '13 at 18:27
  • $\begingroup$ That isn't really a problem. I initialize the table as having zeroes in every element. In order for the Do loop to add to it with each step using the +=, sumlist has to be a list. I suppose I could have it be a list with Null at each element, but I don't think that would make a difference. $\endgroup$
    – user7268
    Jun 4 '13 at 18:28

Another method with new tensor functions introduced in Mathematica 9:

coefficientlist = RandomReal[{0, 10}, {3, 3}];
cubelist = RandomReal[{0, 10}, {3, 5, 5, 5}];

tensor = Fold[
   TensorContract[TensorProduct[#1, cubelist], #2] &,
   {{{1, 3}}, {{1, 5}}}

sum = Outer[tensor[[##, ##]] & @@ {#1, #2, #3} &, Range[5], Range[5], Range[5]]


As I tested with your original scale, both @0x4A4D's method and mine will need lots of memory due to large intermediate arrays. So if you are running short of RAM, like me, Compile your original version might be a better choice.

funcComp = Compile[{{coefficientlist, _Real, 2}, {cubelist, _Real, 4}},
  Module[{sumlist, l},
   sumlist = 0 cubelist[[1]];
   l = Dimensions[coefficientlist][[1]];
   Do[sumlist += 
     coefficientlist[[n, m]] (cubelist[[n]] cubelist[[m]]),
      {n, l}, {m, l}];

Test on a random data:

j = 100;
k = 30;
coefficientlist = RandomReal[{0, 10}, {j, j}];
cubelist = RandomReal[{0, 10}, {j, k, k, k}];

AbsoluteTiming[sum = funcComp[coefficientlist, cubelist];]


sumlist = 0 cubelist[[1]];
 Do[sumlist += 
   coefficientlist[[n, m]] (cubelist[[n]] cubelist[[m]]),
    {n, j}, {m, j}]]


sumlist == sum


  • $\begingroup$ That's great! Originally when you said to use compile, I did something like func=Compile[{{n,_Integer},{m,_Integer}},coefficientlist[[n, m]] (cubelist[[n]] cubelist[[m]])], but didn't see any speedup. $\endgroup$
    – user7268
    Jun 4 '13 at 20:02
  • $\begingroup$ Then, I tried copying from your example, and it took a fraction of a second, and thought "WOW!" before I realized you had the summation running from 1 to 3, LOL. But even doing the full summation, it dropped the time needed from ~1100 seconds to 123 seconds. Thank you very much. $\endgroup$
    – user7268
    Jun 4 '13 at 20:04
  • $\begingroup$ Since both sums have different ranges how can sumlist == sum? You probably should change the ranges and timing of the second sum. $\endgroup$ Jun 4 '13 at 20:08
  • $\begingroup$ @SjoerdC.deVries Sorry for the mistake.. It is indeed equal due to a same initialization and computing process, which could be my excuse of the oversight..(also too late/early to have a clear brain here. :-/ Anyway, please see my edit. $\endgroup$
    – Silvia
    Jun 4 '13 at 20:29
  • $\begingroup$ @user7268 Glad I can help:) Please see my edit on the Compile part. $\endgroup$
    – Silvia
    Jun 4 '13 at 20:32

At the very least, here is a loop-free method:

Total[coefficientlist Outer[Times, cubelist, cubelist, 1], 2]

There surely are more efficient methods, tho.

  • 2
    $\begingroup$ The problem here is that generating the matrix Outer[Times, cubelist, cubelist, 1], which would have dimension {176,176,65,65,65}, takes too much memory and causes the system to crash. Also, does MapThread help this to be more efficient? I would have just tried Total[coefficientlist Outer[Times, cubelist, cubelist, 1],2] since it assumes element by element multiplication anyway. $\endgroup$
    – user7268
    Jun 4 '13 at 19:22
  • 1
    $\begingroup$ Hmm, right; the Hadamard multiplication can be done more directly. I'll edit. On a hunch, what are the results of Developer`PackedArrayQ[cubelist] and Developer`PackedArrayQ[coefficientlist]? $\endgroup$
    – J. M.'s torpor
    Jun 4 '13 at 19:24
  • $\begingroup$ Hah, it does have a name! Thank you for that piece of information, 0x4A4D. $\endgroup$
    – user7268
    Jun 4 '13 at 19:27
  • $\begingroup$ You haven't answered my question yet... did you try executing those two snippets? $\endgroup$
    – J. M.'s torpor
    Jun 5 '13 at 1:01
  • $\begingroup$ Sorry for not answering your question, I did run those and the answer is False. I can't find from the online documentation what PackedArrays are good for. Converting the arrays in question into PackedArrays didn't improve the performance of the solution Silvia posted. $\endgroup$
    – user7268
    Jun 5 '13 at 16:55

Your Answer

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