As an experiment, I tried to write code for simulating many realizations of a homogeneous Poisson point process in a very fast fashion. My idea was to simulate all the realizations in two short steps.
In reality, the findings of this experiment and the contents of this post have little practical value, as computers are so fast at generating Poisson point processes. Still, it was an interesting task, which taught me a couple of things. And I did produce faster code.
MATLAB
I first tried this experiment in MATLAB.
Vectorization
In languages like MATLAB, the trick for speed is to use vectorization, which means applying a single operation to an entire vector or matrix (or array) without using an explicit for-loop. Over the years, the people behind MATLAB has advised to use vectorization instead of for-loops, as for-loops considerably slowed down MATLAB code. (But, as as time goes by, it seems using for-loops in MATLAB doesn’t slow the code down as much as it used to.)
Simulating Poisson point processes is particularly amenable to vectorization, due to the independent nature of the point process. One can simulate the number of points in each realization for all realizations in one step. Then all the points across all realizations can also be positioned in one step. In the two-dimensional case, this results in two one-dimensional arrays (or vectors, in MATLAB parlance) for the \(x\) and \(y\) coordinates. (Of course, in my code, I could have used just one two-dimensional array/vector for the coordinates of the points, but I didn’t.)
After generating the points, the coordinates of the points need to be grouped into the different realizations and stored in appropriate data structures.
Data structures
The problem with storing point processes is that usually each realization has a different number of points, so more sophisticated data structures than regular arrays are needed. For MATLAB, each point process realization can be stored in a data object called a cell array. For a structure array, it’s possible for each element (that is, structure) to be a different size, making them ideal for storing objects like point processes with randomly varying sizes.
In the case of two-dimensional point processes, two cell arrays can be used to store the \(x\) and \(y\) coordinates of all the point process realizations. After randomly positioning all the points, they can be grouped into a cell array, where each cell array element represents a realization of the Poisson point process, by using the inbuilt function MATLAB mat2cell, which converts a matrix (or array) into a cell array.
Alternatively, we could use another MATLAB data object called a structure array. In MATLAB structures have fields, which can be, for example for a point process can be the locations of the points. Given cell arrays of equal size, we can convert them into a single structure array by using the inbuilt MATLAB function struct.
Python
After successfully simulating Poisson point processes in MATLAB, I tried it in Python with the NumPy library.
Vectorization
I basically replicated what I did in MATLB using Python by positioning all the points in a single step. This gives two one-dimensional NumPy arrays for the \(x\) and \(y\) coordinates of all the point process realizations. (Again, I could have instead stored the coordinates as a single two-dimensional array, but I didn’t.)
Perhaps surprisingly, the vectorization step doesn’t speed things up much in Python with the NumPy library. This may be due to the fact that the underlying code is actually written in the C language. That motivated me to see what methods have been implemented for simulating Poisson variables, which is the topic of the next couple posts.
Data structures
In Python, the data structure known as a list is the natural choice. Similarly to cell arrays in MATLAB, two lists can be used for the \(x\) and \(y\) coordinates of all the points. Instead of MATLAB’s function mat2cell, I used the NumPy function numpy.split to create two lists from the two NumPy arrays containing the point coordinates.
Python does not seem to have an immediate equivalent to structure arrays in MATLAB. But in Python one can define a new data structure or class with fields, like a structure. Then one can create a list of those data structures with fields, which are called attribute references in Python.
Code
The code in MATLAB and Python can be found here. For a comparison, I also generated the same number of point process realizations (without using vectorization) by using a trusty for-loop. The code compares the times of taken for implemented the two different approaches, labelled internally as Method A and Method B. There is a some time difference in the MATLAB code, but not much of a difference in the Python case.
I have commented out the sections that create data structures (with fields or attribute references) for storing all the point process realizations, but those sections should also work when uncommented.