Re: Re: Using Fourier and InverseFourier instead of Convolve

*To*: mathgroup at smc.vnet.net*Subject*: [mg86003] Re: [mg85958] Re: Using Fourier and InverseFourier instead of Convolve*From*: Daniel Lichtblau <danl at wolfram.com>*Date*: Sat, 1 Mar 2008 04:43:25 -0500 (EST)*References*: <fpooi0$1cg$1@smc.vnet.net> <fq0u59$k8l$1@smc.vnet.net> <200802280759.CAA19967@smc.vnet.net>

Solomon, Joshua wrote: > First, let me thank all who responded to my query! > Second, no, Daniel, your code works fine for me. > Finally, FWIW here is code that does what I really wanted, which was to > quickly convolve a long list of lists (each having the same length). > Possibly not the most elegant code ever, but it works, as I said, FWIW: > > cv[listOfLists_]:=With[{dims=Dimensions[listOfLists]}, > (dims[[1]]dims[[2]]-dims[[1]]+1)^((dims[[1]]-1)/2) > InverseFourier[Times[##]]&@@( > Fourier[PadRight[#,dims[[1]]dims[[2]]-dims[[1]]+1] > ]&/@listOfLists)] > [...] (Just to note, I am not the Daniel referred to above.) If speed is the interest, and you have many lists, you'll generally do better to use a divide-and-conquer approach. Here is your code, name changed. CauchyMultiply1[listOfLists_] := With[ {dims=Dimensions[listOfLists]}, (dims[[1]]*dims[[2]]-dims[[1]]+1)^((dims[[1]]-1)/2)* InverseFourier[Times[##]]& @@( Fourier[PadRight[#,dims[[1]]*dims[[2]]-dims[[1]]+1] ]&/@listOfLists)] Here is similar but using ListConvolve iteratively. CauchyMultiply2[listOfLists_] := Fold[ ListConvolve[#1, #2, {1,-1}, 0]&, {1}, listOfLists] Here is a divide-and-conquer using ListConvolve. CauchyMultiply3[ll_] := With[{len=Length[ll]}, If [len==1, ll[[1]], If [len==2, ListConvolve[ll[[1]], ll[[2]], {1,-1}, 0], ListConvolve[ CauchyMultiply3[Take[ll,Floor[len/2]]], CauchyMultiply3[Drop[ll,Floor[len/2]]], {1,-1}, 0]]]] I will show an example on lists of integers. This makes the second method quite slow (because it is exact). For a list of machine floats it tends to be faster than the first (yours) by a factor of 2 or so. In[84]:= lists = RandomInteger[{-10,10}, {100,100}]; In[85]:= Timing[cm1 = CauchyMultiply1[lists];] Out[85]= {0.568035, Null} In[86]:= Timing[cm2 = CauchyMultiply2[lists];] Out[86]= {19.4332, Null} In[87]:= Timing[cm3 = CauchyMultiply3[lists];] Out[87]= {0.844053, Null} At this point we'll check that either relative or absolute error is small in eavery position. In[90]:= InputForm[Max[Abs[(cm1-cm3)/Max[1,Sqrt[cm1^2+cm3^2]]]]] Out[90]//InputForm= 3.645529318632443*^-15 Also note that cm2 and cm3 exactly agree (as they should, since they both use integer arithmetic). In[91]:= Max[Abs[cm2-cm3]] Out[91]= 0 Thus far I've shown no reason to use CauchyMultiply3, if your inputs are approximate. We'll see now that it is quite fast. In[92]:= Timing[cm4 = CauchyMultiply3[N[lists]];] Out[92]= {0.024003, Null} The error is still within reason. In[96]:= InputForm[Max[Abs[(cm4-cm3)/Max[1,Sqrt[cm4^2+cm3^2]]]]] Out[96]//InputForm= 1.039317649229334*^-14 But it was around 20x faster than CauchyMultiply1. I will add that the memory hoofprint of CauchyMultiply3 is also typically smaller. Below is from a different session. I use exact integer arithmetic for this example, and in interest of full disclosure I will mention that CauchyMultiply3 does not handle it correctly at machine precision (I do not know why, but probably a limitation somwhere in Fourier; I plan to file a report on this). In[6]:= lists = RandomInteger[{-10,10}, {1000,10}]; In[7]:= Timing[cm3 = CauchyMultiply3[lists];] Out[7]= {5.27, Null} In[8]:= MemoryInUse[] Out[8]= 9809512 In[9]:= MaxMemoryUsed[] Out[9]= 40200072 In[10]:= Timing[cm1 = CauchyMultiply1[lists];] Out[10]= {23.56, Null} In[11]:= MemoryInUse[] Out[11]= 10427448 In[12]:= MaxMemoryUsed[] Out[12]= 656384928 So we got an exact result around 4.5 times faster than the approximate one, and max memory required was 40 vs 656 meg. Observe that the memory needed for the results goes in the opposite direction, because the exact result uses bignum integers instead of machine floats. Daniel Lichtblau Wolfram Research