Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

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
2  
@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
1  
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

2 Answers 2

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] &,
   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

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
2  
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
1  
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

 
discard

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