3

Is there a "Mathematica Way", like Map or Apply to compute the following double sum?

$\sum_{i=1}^{N_1}\sum_{j=1}^{N_2} m_i n_j \, f(\tau_{i} \gamma_{j})$

I have already stored the lists $m,n,\tau,\gamma$.

$m$ and $\tau$ are of length $N_1$ and $n$ and $\gamma$ of length $N_2$. So my approach was to use something like

Outer[f, τ, γ]

for the function map, and then 

Total[Total[Outer[f, τ, γ]]]

but how do Ι multiply mand n into the sum?

EDIT: Undo Edit, figured out the bug in my code

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
rainer
  • 401
  • 4
  • 12
  • I guess the time to beat is Sum[m[[i]] n[[j]] f[tau[[i]], gamma[[j]]], {i, 1, n1}, {j, 1, n2}] – ssch Feb 05 '13 at 16:34
  • If f is simple enough then Compile will likely give the best possible performance. But it's not "Mathematica-like". – Szabolcs Feb 05 '13 at 16:41
  • You could likewise do an Outer[Times,m,n] and then Dot that with your first Outer result. I don't know how that will be speedwise though. – Daniel Lichtblau Feb 05 '13 at 16:45
  • Does Table/ParallelTable count as Mathematica-like? – RunnyKine Feb 05 '13 at 16:50
  • ;) Thank you for your answers. For me Mathematica-like is to avoid for or while loops and to use fast built-in functions – rainer Feb 05 '13 at 17:06
  • this seems a miniscule bit faster..: Flatten[ Table[ m[[i]] n[[j]] , {i, n1}, {j, n2}] ].Flatten[Table[f[t[[i]],z[[j]]], {i, n1}, {j, n2}]]. – george2079 Feb 05 '13 at 18:01
  • 1
    Actually I missed the need for Flatten. This variant works. AbsoluteTiming[Flatten[Outer[Times, m, n]].Flatten[ Outer[f, t, gamma]]] – Daniel Lichtblau Feb 05 '13 at 20:36

3 Answers3

9

If your elements are in lists the fastest way is to use array operations. In the present case of an outer product one index, let's say "i", will not be expanded, on the other you want to thread.

To operate on a list the function needs the attribute Listable. The Times function, as many other internal ones, is already listable, that is

Times[{1,2,3},x] = {x, 2x, 3x}

You want to do the same with your function "f", e.g.

ClearAll[f]
SetAttributes[f, Listable]

f[{1,2,3},x] = {f[1, x], f[2, x], f[3, x]}

For the threading part you can use MapThread to "slide" through the list, e.g.

MapThread[#1 f[#2] &, {{a,b,c},{1,2,3}}] = {a f[1], b f[2], c f[3]}

Finally your code will look like this

ClearAll[f]
SetAttributes[f, Listable]
f[x_, y_] := x y
{m, n, t, g} = RandomReal[1, {4, 1000}];

Total@Total@MapThread[m #1 f[t #2] &, {n, g}]

A note about speed: usually threaded operation are much faster than element ones

In[27]:= 
Total@Total@MapThread[m #1 f[t, #2] &, {n, g}] // AbsoluteTiming
Total@Total[Outer[Times, m, n] Outer[f[#1, #2] &, t, g]] // AbsoluteTiming
Sum[m[[i]] n[[j]] f[t[[i]], g[[j]]], {i, 1, 1000}, {j, 1, 1000}] // AbsoluteTiming
ParallelSum[m[[i]] n[[j]] f[t[[i]], g[[j]]], {i, 1, 1000}, {j, 1, 1000}] // AbsoluteTiming

Out[27]= {0.797645, 62806.3}

Out[28]= {1.554886, 62806.3}

Out[29]= {2.933898, 62806.3}

Out[30]= {0.828527, 62806.3}

Thus always try to use them when possible.

Batracos
  • 2,105
  • 14
  • 14
1

How about using Table or ParallelTable like this:

Total[ParallelTable[m[[i]] n[[j]] f[τ[[i]] γ[[j]]], {i, 1, N1}, {j, 1, N2}]]
RunnyKine
  • 33,088
  • 3
  • 109
  • 176
1

You almost had it.You were right to use Outer[]; what you missed was to interpret the sum as an appropriate matrix multiplication. Observe:

Table[Subscript[m, i], {i, 5}].Outer[f, Table[Subscript[τ, i], {i, 5}],
                                        Table[Subscript[γ, j], {j, 4}]].
                               Table[Subscript[n, j], {j, 4}] == 
Sum[f[Subscript[τ, i], Subscript[γ, j]] Subscript[m, i] Subscript[n, j],
    {i, 5}, {j, 4}] // FullSimplify
   True
J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574