Faster code for specific problem

Question

Sorry, but I need the following code much faster and I'm running out of ideas. I do appreciate any help. I post two of my fastest ideas (tab and mat), although they are not really fast. Thanks in advance.

Clear[bNum,nNum,weights,newMat,sortMat,tab];

bNum=20;
nNum=100;

weights=Table[RandomReal[{-10,10}],{b,1,bNum},{n,1,nNum}];
newMat=Table[RandomReal[{-10,10}],{b,1,bNum},{n,1,nNum}];
sortMat=Table[i,{i, -10, 10, 0.1}];

tab=Table[
Sum[weights[[k,i]]*If[newMat[[k,i]] <= sortMat[[j]],1.,0.],{i,1,nNum}],{k,1,bNum}, 
{j,1,Length[sortMat]}];//AbsoluteTiming

mat=Block[{k=1},
Reap[Do[
Sow[Table[Sum[weights[[k,i]]*If[newMat[[k,i]] <= sortMat[[j]],1.,0.],{i,1,nNum}],{j,1,Length[sortMat]}]];
k++,{bNum}]]][[2,1]];//AbsoluteTiming

Answer 1

As I have advocated in another thread, usually the best way (by far) to speed up a calculation is first to improve the algorithm. Then--if it remains necessary--you can focus on platform-specific methods to speed up execution.

Here, newMat and weights are viewed as lists of vectors to be processed in parallel and (according to a comment to the question) sortMat is really a sequence $(y_0, y_0 + dy, y_0 + 2 dy, \ldots, y_0 + (n-1) dy)$. For each $y$ in this sequence, the algorithm seeks the sum of weights associated with elements $x$ of newMat less than or equal to $y$.

[To find those elements, we change the units of measurement of x and sortMat so that the new sortMat is the integral sequence $(1, 2, \ldots, n)$ and we sort x (and sort the weights in parallel with x to keep their association intact). After these preliminaries, the integer part of any entry in x tells us precisely how many elements of sortMat are less than it. Thus--work a few simple examples to check--the differences in this sorted, re-expressed x tell us how many times each cumulative weight ought to appear in the output. The larger sortMat is relative to x, the more computation this approach saves.]

An efficient approach would capitalize on a preliminary sorting of newMat (which requires sorting the corresponding vector weights in parallel). Do this with Ordering. This reduces the amount of computation to the length of each row of newMat rather than the length of sortMat and it reduces a triple nested loop to a double nested loop (which is the minimum possible, given that the output is a rank 2 array).

ClearAll[process];
process[x_, w_, {y0_, dy_, ny_}] := Module[{x0, w0, v, i, f},
   x0 = Min[ny, Max[0, Ceiling[(# - y0)/dy]]] & /@ x;
   x0 = {0}~Join~Sort[x0]~Join~{ny};
   w0 = {0}~Join~Accumulate[w[[Ordering[x]]]];
   f := Function[{v, i}, Flatten[ConstantArray[#[[1]], #[[2]]] & /@ ({v, i}\[Transpose])]];
   f[w0, Differences[x0]]
   ];

Let's check the timing and the correctness of this solution. Because it's much faster, let's make newMat have ten times as many rows as before:

bNum = 200;
nNum = 100;
weights = Table[RandomReal[{-10, 10}], {b, 1, bNum}, {n, 1, nNum}];
newMat = Table[RandomReal[{-10, 10}], {b, 1, bNum}, {n, 1, nNum}];
sortMat = Table[i, {i, -10, 10, 0.1}];

Now:

tab = Table[Sum[weights[[k, i]]*
      If[newMat[[k, i]] <= sortMat[[j]], 1., 0.], {i, 1, nNum}], 
          {k, 1, bNum}, {j, 1, Length[sortMat]}]; // AbsoluteTiming

{8.2744732, Null}

and (describing sortMat with the list $(y_0, dy, n)$ = $(-10, 0.1, 201)$)

tab0 = MapThread[process[#1, #2, {-10, 0.1, 201}] &, {newMat, weights}]; // AbsoluteTiming

{0.0920053, Null}

It's 89 times faster, yet produces the same result to within floating point error:

Plus @@ (Abs[Chop[#]] & /@ (Flatten[tab] - Flatten[tab0])) == 0.0

True

Now, if you wish to speed up the improved algorithm with a compilation, go ahead.

---

Edit

The comments and a further chat session (linked in the comments) uncovered some implicit assumptions in the original solution and led to an improvement that more closely matches the original. The two algorithms are not completely identical when applied to machine-precision numbers, though, due to differences in how the values in sortMat are computed and compared to the values in newMat. Testing indicates they do agree except when elements of newMat appear to equal elements of sortMat: an equality test might or might not return True in such cases.

Answer 2

Updated with whuber's fix, and added process3

Murta showed Compile and whuber showed a great algorithm improvement.
Let me show some syntax tweaking.

whuber's process function can be written more efficiently and concisely like this:

process2[x_, w_, {y0_, dy_, ny_}] :=
  Module[{x0, w0},
    x0 = Ceiling[(x - y0)/dy] ~Clip~ {0, ny};
    x0 = Join[{0}, Sort @ x0, {ny}];
    w0 = {0} ~Join~ Accumulate[ w[[Ordering @ x]] ];
    Join @@ MapThread[ConstantArray, {w0, Differences@x0}]
  ]

Your code producing the random tables can be made much more efficient, as there is a specific syntax within RandomReal (and RandomInteger, etc.) for this.

bNum = 200;
nNum = 100;
weights = RandomReal[{-10, 10}, {bNum, nNum}];
newMat  = RandomReal[{-10, 10}, {bNum, nNum}];

Timing comparison with process:

MapThread[process[#1, #2, {-10, 0.1, 201}] &, {newMat, weights}]; //Timing //First

MapThread[process2[#1, #2, {-10, 0.1, 201}] &, {newMat, weights}]; //Timing //First

0.0842

0.01744

(My code originally showed about a doubling of performance over whuber's. It's now nearly 5X due to the poor performance of Min[ny, Max[0, Ceiling[(# - y0)/dy]]] & /@ x in the fix, which I handled with Clip.)

---

I decided to implement a form of the function that operates on newMat and weights directly, that is moving the outer MapThread inside the function. This allows certain optimizations at the cost of additional memory consumption.

process3[x_, w_, {y0_, dy_, ny_}] :=
  Module[{xm, wm},
    xm = Ceiling[(x - y0)/dy] ~Clip~ {0, ny};
    xm = ArrayFlatten @ {{0, Sort /@ xm, ny}};
    wm = ArrayFlatten @ {{0, MapThread[Accumulate @ #2[[Ordering@#]] &, {x, w}]}};
    Join @@ MapThread[ConstantArray, {#2, Differences@#}] & ~MapThread~ {xm, wm}
  ]

process3[newMat, weights, {-10, 0.1, 201}] //Timing //First

0.01496

Answer 3

Here is one demonstration on how to optimize your code using compile. I see that you are not using the functional programin mindset yet, this could be improved in your code too, but here I will focus just on compile.

Here we create your original matrices:

Clear[bNum,nNum,weights,newMat,sortMat,tab];

bNum=20;
nNum=100;

weights = RandomReal[{-10, 10}, {bNum, nNum}];
newMat = RandomReal[{-10, 10}, {bNum, nNum}];
sortMat = Range[-10, 10, 0.1];

Now I created 3 exemples of your tab computation:

tab[weights_,newMat_,sortMat_]:=Module[{bNum,nNum},
    {bNum,nNum}=Dimensions[weights];
    Table[
            Sum[
                weights[[k,i]]*If[newMat[[k,i]]<=sortMat[[j]],1.,0.]
                ,{i,1,nNum}
                ]
            ,{k,1,bNum},{j,1,Length[sortMat]}
    ]
]

tabCompiled=Compile[{{weights,_Real,2},{newMat,_Real,2},{sortMat,_Real,1}},

    Module[{bNum,nNum},
    {bNum,nNum}=Dimensions[weights];
    Table[
            Sum[
                weights[[k,i]]*If[newMat[[k,i]]<=sortMat[[j]],1.,0.]
                ,{i,1,nNum}
                ]
            ,{k,1,bNum},{j,1,Length[sortMat]}
        ]
    ]
];

tabCompiledC=Compile[{{weights,_Real,2},{newMat,_Real,2},{sortMat,_Real,1}},

    Module[{bNum,nNum},
    {bNum,nNum}=Dimensions[weights];
    Table[
            Sum[
                weights[[k,i]]*If[newMat[[k,i]]<=sortMat[[j]],1.,0.]
                ,{i,1,nNum}
                ]
            ,{k,1,bNum},{j,1,Length[sortMat]}
        ]
    ]
    ,CompilationTarget->"C"
];

The first without compile, the second compiled, the third compiled in C. Compare the diferences:

tab[weights, newMat, sortMat]; // AbsoluteTiming
tabCompiled[weights, newMat, sortMat]; // AbsoluteTiming
tabCompiledC[weights, newMat, sortMat]; // AbsoluteTiming
{0.833, Null}
{0.046, Null}
{0.016, Null}

I hope that it could help you.

Answer 4

You can get a decent speed up with a little tweak of the method suggested in Daniel Lichtblau's comment, by transposing newMat and weights so that Total can be used more efficiently:

tab2 = Module[{nt = Transpose[newMat], wt = Transpose[weights]}, 
  Transpose[Total[UnitStep[# - nt] wt] & /@ sortMat]];

With the array sizes in the question this gives me a speedup factor of about 120.

The code can be dropped directly into Compile for another factor of 2:

c = Compile[{{newMat, _Real, 2}, {weights, _Real, 2}, {sortMat, _Real, 1}},
   Module[{nt = Transpose[newMat], wt = Transpose[weights]}, 
    Transpose[Total[UnitStep[# - nt] wt] & /@ sortMat]]];

tab3 = c[newMat, weights, sortMat];