MathGroup Archive 2008

[Date Index] [Thread Index] [Author Index]

Search the Archive

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




  • Prev by Date: Mathematica 6 obtains imaginary eigenvalues for a Hermitian matrix
  • Next by Date: Re: Finding a continuous solution of a cubic
  • Previous by thread: Re: Mathematica 6 obtains imaginary eigenvalues for a Hermitian matrix
  • Next by thread: Re: Finding a continuous solution of a cubic