Can I do this faster?
in [Mathematica]
|
Prev: plot question
Next: Faster alternative to AppendTo? --
From: Andreas on 3 Sep 2009 05:30 I start with two 2 dimensional lists. The have equal depths, but may or may not have equal lengths. Possible data looks like this: list1 = RandomReal[{0, 1}, {1500, 3}]; list2 = RandomReal[ExponentialDistribution[2], {1000, 3}]; The following produces the results I need, but takes longer to run than I'd like: Rest[FoldList[Times, 1, Transpose[(t = #; 1 + Inner[Times, t, # - 1, Plus] & /@ list1) & /@ list2]]] Any suggestions rewriting this so I can make it run faster? Thanks.
From: Szabolcs Horvát on 3 Sep 2009 19:54 On 2009.09.03. 11:30, Andreas wrote: > Rest[FoldList[Times, 1, Transpose[(t = #; 1 + Inner[Times, t, # - 1, > Plus]& /@ list1)& /@ > list2]]] First I'd like to say that it's often much easier (at least for me!) to come up with a solution if you explain what you want to do in plain English instead of just providing a program to speed up. Now I have to convert the program to a form that my brain can handle, then convert it back to program code... Why not avoid the first step if possible? ;-) So, a few things we might notice about this implementation: 1. Inner is used with Times and Plus, so why not replace it with Dot? 2. That nested function (with the assignment to the global t) looks discomforting. I'm not sure how Mathematica's compiler can handle that (Map auto-compiles the function when working on large lists). So let's try to get rid of that also. These might not be the main reson for the slowdown. I am just guessing---predicting Mathematica's performance can be difficult. So, rewrite the inner part of the program first. Instead of Inner we can use Dot, instead of the nested function we can use Outer: list1 = RandomReal[{0, 1}, {1500, 3}]; list2 = RandomReal[ExponentialDistribution[2], {1000, 3}]; In[3]:= Timing[ x = Transpose[(t = #; 1 + Inner[Times, t, # - 1, Plus] & /@ list1) & /@ list2];] Out[3]= {21.656, Null} In[4]:= Timing[y = Outer[1 + #2.(#1 - 1) &, list1, list2, 1];] Out[4]= {10.625, Null} That's a 2x speedup. Are the results equivalent? In[5]:= x == y Out[5]= False We got False, but that's only because of numerical errors (the operations are performed in a different order): In[6]:= Max(a)Abs[x - y] Out[6]= 8.88178*10^-16 So the result is correct. What else can we do to speed things up? Notice that it is not necessary to subtract/add 1 in the inner function 1 + #2.(#1 - 1) &. This can be done on the input and output instead. So we can get rid of custom functions and use the built-in Dot only: In[7]:= Timing[z = 1 + Outer[Dot, list1 - 1, list2, 1];] Out[7]= {1.25, Null} In[8]:= y == z Out[8]= True Now that's a 17x speedup compared to the implementation we started with. Trying to simplify things will often pay off because it will be easier to see how to rewrite the program to use built-in functions and packed arrays as much as possible. It's also much easier to see what the program does. The FoldList part takes an additional 3 seconds on my machine. I can't help with speeding that up unfortunately. I hope this helps, Szabolcs
From: Daniel Lichtblau on 3 Sep 2009 19:55 Andreas wrote: > I start with two 2 dimensional lists. The have equal depths, but may > or may not have equal lengths. Possible data looks like this: > > list1 = RandomReal[{0, 1}, {1500, 3}]; > list2 = RandomReal[ExponentialDistribution[2], {1000, 3}]; > > The following produces the results I need, but takes longer to run > than I'd like: > > Rest[FoldList[Times, 1, Transpose[(t = #; 1 + Inner[Times, t, # - 1, > Plus] & /@ list1) & /@ > list2]]] > > Any suggestions rewriting this so I can make it run faster? > > Thanks. Could do: Rest[FoldList[Times, 1, (list1-1).Transpose[(list2)]+1]] One bottleneck is that you will run into software arithmetic because numbers become too small for machine double representation. If you don't mind having them clipped, I think the code below might do what you want. foldTimesC = Compile[{{mat,_Real,2}}, Module[ {current=ConstantArray[1.,Length[mat[[1]]]]}, Table[current = current*mat[[j]], {j,Length[mat]}]]] foldTimesC[(list1-1).Transpose[(list2)]+1]; Daniel Lichtblau Wolfram Research
From: Ray Koopman on 3 Sep 2009 19:56 Rest(a)FoldList[Times,1,1+(list1-1).Transpose(a)list2] On Sep 3, 2:30 am, Andreas <aa...(a)ix.netcom.com> wrote: > I start with two 2 dimensional lists. The have equal depths, but may > or may not have equal lengths. Possible data looks like this: > > list1 = RandomReal[{0, 1}, {1500, 3}]; > list2 = RandomReal[ExponentialDistribution[2], {1000, 3}]; > > The following produces the results I need, but takes longer to run > than I'd like: > > Rest[FoldList[Times, 1, Transpose[(t = #; 1 + Inner[Times, t, # - 1, > Plus] & /@ list1) & /@ list2]]] > > Any suggestions rewriting this so I can make it run faster? > > Thanks.
From: pfalloon on 4 Sep 2009 03:12
On Sep 3, 7:30 pm, Andreas <aa...(a)ix.netcom.com> wrote: > I start with two 2 dimensional lists. The have equal depths, but may > or may not have equal lengths. Possible data looks like this: > > list1 = RandomReal[{0, 1}, {1500, 3}]; > list2 = RandomReal[ExponentialDistribution[2], {1000, 3}]; > > The following produces the results I need, but takes longer to run > than I'd like: > > Rest[FoldList[Times, 1, Transpose[(t = #; 1 + Inner[Times, t, # - 1, > Plus] & /@ list1) & /@ > list2]]] > > Any suggestions rewriting this so I can make it run faster? > > Thanks. I think the main reason for the slowness is that the Inner function is much less efficient for large numerical matrices than the standard Dot function. When you re-express the function using Dot, it is much faster: (* original definition *) f1[list1_,list2_] := Module[{t}, Rest[FoldList[Times, 1, Transpose[(t = #; 1 + Inner[Times, t, # - 1, Plus] & /@ list1) & /@ list2]]]] (* first attempt: incorporate the addition & subtraction into the Inner function *) f2[list1_,list2_] := Inner[Times[#1-1,#2] &, list1, Transpose(a)list2, Plus[1, ##] &] // FoldList[Times, 1, #] & // Rest (* best attempt: apply the addition & subtraction to the matrices as a whole, so we can use the more-efficient Dot function *) f3[list1_,list2_] := Rest(a)FoldList[Times, 1, 1+(list1-1).Transpose [list2]] Here is an example to show the speed-up: In[62]:= (* sample matricies *) list1 = RandomReal[{0, 1}, {1000, 3}]; list2 = RandomReal[ExponentialDistribution[2], {1500, 3}]; In[64]:= (* timings *) res1=f1[list1, list2]; // AbsoluteTiming res2=f2[list1, list2]; // AbsoluteTiming res3=f3[list1, list2]; // AbsoluteTiming Out[64]= {14.0311602,Null} Out[65]= {9.7811874,Null} Out[66]= {1.7187390,Null} However, in testing this I noticed that the function appears to be highly unstable numerically -- so you may want to investigate that further: (* test correctness of results: appears to be numerically unstable! *) Max /@ Abs[res1 - res2] Max /@ Abs[res1 - res3] Hope this helps! Cheers, Peter. |