# 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? ;} –  David Carraher 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

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[5], Range[5], Range[5]]


# 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[[1]];
l = Dimensions[coefficientlist][[1]];
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[[1]];
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. –  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.

-
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