Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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.

share|improve this question
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... – Rod Jun 4 '13 at 18:19
@Rod Lm The sum is zero even if the matrix is really, really large? ;} – DavidC Jun 4 '13 at 18:26
@RodLm I think sumlist is a container to record the sum, sumlist=Table[0,{65},{65},{65}] is merely an initialization. – Silvia Jun 4 '13 at 18:27
Another question: if cubelist is three-dimensional, how can it have Dimension {176,65,65,65}? – Rod Jun 4 '13 at 18:27
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. – user7268 Jun 4 '13 at 18:28
up vote 3 down vote accepted

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


share|improve this answer
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. – user7268 Jun 4 '13 at 20:02
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. – user7268 Jun 4 '13 at 20:04
Since both sums have different ranges how can sumlist == sum? You probably should change the ranges and timing of the second sum. – Sjoerd C. de Vries Jun 4 '13 at 20:08
@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. – Silvia Jun 4 '13 at 20:29
@user7268 Glad I can help:) Please see my edit on the Compile part. – 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.

share|improve this answer
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. – user7268 Jun 4 '13 at 19:22
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]? – J. M. Jun 4 '13 at 19:24
Hah, it does have a name! Thank you for that piece of information, 0x4A4D. – user7268 Jun 4 '13 at 19:27
You haven't answered my question yet... did you try executing those two snippets? – J. M. Jun 5 '13 at 1:01
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. – user7268 Jun 5 '13 at 16:55

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.