# Achieve performance that is competitive with numpy [duplicate]

What are the best practices to ameliorate time performance of Mathematica? As a practical example how we can decrease timing below?

Timing[Length[Select[Table[Random[], {10000000}], # < 0.5 &]]]
(*{14.820095, 5001699}*)

Table[Timing[
Length[Select[Table[Random[], {1000000}], # < 0.5 &]]][[1]], {10}]

(*{1.435209, 1.466409, 1.482009, 1.482009, 1.497610, 1.482009, \
1.466409, 1.466409, 1.482009, 1.528810}*)


For instance in less than three seconds Python (NumPy) can execute a whole for instruction, whereas Mathematica needs almost 15 sec to do only one evaluation.

import time

import numpy as np
from numpy.random import random as rng
t = time.process_time()

for i in range(10):
samples=rng(10000000)
flips=(samples<0.5)
np.sum(flips)

elapsed_time = time.process_time() - t
print(elapsed_time)

---->3.4320220000000035

• Random[] has been deprecated, perhaps you mean RandomReal[1, {10000000}]? – Feyre Feb 26 '17 at 11:34
• Thanks for the comment. I forgot it (and Mathematica did not complain:-)!). In any case there is not any considerable improvement. – Dimitris Feb 26 '17 at 11:38
• I am concurring with Szabolcs in marking this question as a duplicate. If you feel that your underlying question is not answered there please edit it to clarify how it is different and provide a different example. Also be aware that the nature of Mathematica is such that it often will not have the time performance on execution that other platforms have; however one often more than makes up for this in the "time performance" of developing the program itself. – Mr.Wizard Feb 26 '17 at 12:16
• Ok, I will accept the duplicate nature of my question. I learnt a lot of new things today and I expand my little knowledge of vectorisation. But before closing the question and returning to the very first phrase of the post: "What are the best practices to improve time performance of Mathematica?" in general (and not only in this specific example)? – Dimitris Feb 26 '17 at 12:27
• @dimitris Pushing things into pure numeric operations is a basic guideline; Szabolcs's package is a great example of this. If that were the focus of your question it would be a duplicate of (29349). Also for less experienced users I would recommend my self-Q&A: (7924) – Mr.Wizard Feb 26 '17 at 13:23

Vectorization is one of the most effective ways to increase performance. The numpy code you show is fast because it uses vectorization.

Vectorization means working with entire arrays instead of element by element: using array arithmetic, array comparisons, etc. Array operations can be implemented very efficiently using SIMD processing, and are also straightforwardly parallelizable. In fact Mathematica uses multiple CPU cores for vector arithmetic. This happens to be important for correct benchmarking as well because Timing adds up the time spent by individual CPU cores while AbsoluteTiming measures how much time has actually elapsed. Thus below I use only AbsoluteTiming. Timing would give longer and inaccurate times (try it!).

The first thing we can do is replace Table[Random[], n] by RandomReal[1, n]. In the spirit of vectorization, generate the whole array in one go instead of element by element:

RandomReal[1, 10000000]; // AbsoluteTiming
(* {0.083058, Null} *)

Table[Random[], 10000000]; // AbsoluteTiming
(* {0.312061, Null} *)


Then instead of Select, use the techniques I described here:

Mathematica does not have built-in vector comparison, like samples < 0.5 in numpy. But it is always possible to express these operations in terms of simple arithmetic, as described in the thread linked above. Unfortunately, this often results in expressions which are hard to decipher. To make it easier to use these techniques, I wrote a small package called BoolEval that will translate expressions written in terms of relational operators (like <, >, ==, etc.) into vector arithmetic.

This is how we can apply the package to your example:

Select[Table[Random[], {10000000}], # < 0.5 &]; // AbsoluteTiming
(* {5.02988, Null} *)

<< BoolEval

AbsoluteTiming[
arr = RandomReal[1, 10000000];
BoolPick[arr, arr < 0.5];
]
(* {0.310938, Null} *)


Behind the scenes, BoolEval translates this into:

Pick[arr, 1 - UnitStep[arr - 0.5], 1]; // AbsoluteTiming
(* {0.237126, Null} *)


We can see this by using BoolEval with a symbolic expression:

BoolEval[a < 0.5]
(* 1 - UnitStep[-0.5 + a] *)


The package also has other useful functions such as BoolCount for counting how many elements satisfy the condition:

BoolCount[arr < 0.5] // AbsoluteTiming
(* {0.07805, 4998908} *)


This translates to Total[1 - UnitStep[arr - 0.5]] and is precisely equivalent to the numpy code you showed. As a comparison, your numpy code runs in ~0.14 seconds on my machine, which is nearly twice slower. I was using an MKL-enabled numpy (from Anaconda), so the comparison with Mathematica is not unfair.

• Thank you very much. A lot of homework today:-)! – Dimitris Feb 26 '17 at 12:30
• @dimitris I have Python 2.7, so I used time.clock() instead of time.process_time(). Do you know if that can be a problem? Is it still accurate? (I'm thinking of how Timing is not accurate here, but AbsoluteTiming is.) – Szabolcs Feb 26 '17 at 12:57
• I don't know. I learnt about the time.process_time() today through googling. See this question and the replies given therein stackoverflow.com/questions/7370801/… – Dimitris Feb 26 '17 at 13:03

As noted by Szabolcs, there is an inefficiency in generating the table of real elements one at a time.

v = RandomReal[{0, 1}, 1000000];
Table[Timing[Length[Select[v, # < 0.5 &]]][[1]], {10}]
(* {0.384, 0.38, 0.376, 0.38, 0.38, 0.38, 0.504, 0.48, 0.384, 0.38} *)


However, a more important time saving can be made by compiling the selection function

sel = Compile[{{x, _Real, 1}}, Select[x, # < 0.5 &]];
Table[Timing[Length[sel[v]]][[1]], {10}]
(* {0.072, 0.068, 0.068, 0.068, 0.072, 0.068, 0.064, 0.068, 0.068, 0.068} *)
`
• Compilation is a very good way to speed up such operations (+1), but in this particular case vectorization outperforms it. For an array of length 10,000,000, I get timings of 0.2 s from the vectorized version and 0.75 from the compiled version. I am making this comment only because you are using a 10x shorter array than the OP, thus you are getting 10x shorter timings (0.075 instead of 0.75). This can be a bit misleading when people read the answer. – Szabolcs Feb 26 '17 at 11:52
• @Szabolcs, good point. I cut and pasted the 2nd expression from the original question, which uses the smaller array. – mikado Feb 26 '17 at 11:54
• You are right, I did not notice that. – Szabolcs Feb 26 '17 at 11:57