# Most efficient way to sum over matrix products?

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=Table[0,{65},{65},{65}];
Do[
If[coefficientlist[[n,m]]!=0,
sumlist += coefficientlist[[n,m]] (cubelist[[n]]cubelist[[m]]) ];
,{n,176},{m,176}];


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.

• 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? ;} 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. 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

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] &,
coefficientlist,
{{{1, 3}}, {{1, 5}}}
];

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


# Edit

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[];
l = Dimensions[coefficientlist][];
Do[sumlist +=
coefficientlist[[n, m]] (cubelist[[n]] cubelist[[m]]),
{n, l}, {m, l}];
sumlist
]]


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];]


{4.368250,Null}

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


{7.961455,Null}

sumlist == sum


True

• 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. 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. Jun 4 '13 at 20:29
• @user7268 Glad I can help:) Please see my edit on the Compile part. 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.

• 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 DeveloperPackedArrayQ[cubelist] and DeveloperPackedArrayQ[coefficientlist]? 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? 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