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