Faster code for specific problem

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
• Rephrase your question to be an actual question, format your code (press the golden question mark in the top right of the input text field to read about formatting, it leads to mathematica.stackexchange.com/editing-help ), and explain what you are trying to do rather then just post a code snippet with the request of "make it faster". – jVincent Oct 5 '12 at 14:10
• To restate jVincent's plea: please state your actual problem, format your code, and then we can be more helpful... – J. M. will be back soon Oct 5 '12 at 14:23
• Hi Chris, and welcome to the Mathematica StackExchange. Sorry if some of the first responses you received seem overly critical. We sometimes forget that new folks need some time to understand just how the site works. Don't forget to upvote good answers (and other people's questions) using the triangle above the number next to the post, and use the checkmark to "accept" the answer to your question that you think best answers it. – DavidC Oct 5 '12 at 14:48
• This is faster: Timing[tab2 = Transpose[Table[ Map[Total,weights*UnitStep[sortMat[[j]]-newMat]], {j,Length[sortMat]}]];] In general you want to use Mathematica's list-based operations for this sort of thing. – Daniel Lichtblau Oct 5 '12 at 16:46
• @whuber: yes the list contains always evenly spaced numbers. – chris Oct 6 '12 at 5:52

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

• Might as well do FoldList[Plus, 0, w[[Ordering@x]]] instead of {0} ~Join~ Accumulate[w[[Ordering@x]]]... :) – J. M. will be back soon Oct 8 '12 at 3:29
• @J.M. prior experience suggests that will be slower, but I'll test it now. Yup, slower, at least on version 7. – Mr.Wizard Oct 8 '12 at 3:30
• BTW: you don't need to generate sortMat anymore. :) – J. M. will be back soon Oct 8 '12 at 3:50
• @Mr.Wizard: Thank you very much for another great improvement. – chris Oct 8 '12 at 8:28
• Thank you: collectively, we're now up to about a factor of 150 improvement over the original. Perhaps the last challenge remaining would be to compile this code, thereby combining both ideas. Compile doesn't like the last step (of joining constant arrays of varying lengths). If there is a way to do this efficiently and get it compiled, possibly the timing can be improved further (but probably not dramatically). – whuber Oct 8 '12 at 22:35

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[#[], #[]] & /@ ({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.

• I'd replace If[newMat[[k, i]] <= sortMat[[j]], 1., 0.] with Boole[newMat[[k, i]] <= sortMat[[j]]]. – J. M. will be back soon Oct 8 '12 at 3:30
• @J.M. You do see that's in the non-optimized part, right? – Mr.Wizard Oct 8 '12 at 3:33
• @Mr. Wizard: Ah, right. After looking at the method in more detail, I think a judicious use of Dot[] might lead to a better route... – J. M. will be back soon Oct 8 '12 at 3:36
• @J.M. I'd like to see another good use of Dot. Please. – Mr.Wizard Oct 8 '12 at 3:39
• @whuber: thanks a lot for this impressive speed improvement. This is exactly what I was looking for. Thanks to all the others as well. – chris Oct 8 '12 at 8:24

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.

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];