Vectorization and arrays: Period-doubling bifurcations and the Mandelbrot Set

In this notebook, we’ll cover the basics of vectorization and arrays, including how to use them to compute the period-doubling bifurcations of the logistic map, and to generate the Mandelbrot set.

External libraries¶

A major benefit of the Python ecosystem is the ability to easily import external libraries. This is done with the import command. For example, the following imports the numpy library, which contains many useful functions for scientific computing.

import numpy as np

Common libraries you might use in the natural sciences include

• matplotlib for general plotting
• numpy for linear algebra and lots of MATLAB functions
• scipy for common scientific algorithms
• scikit-learn for machine learning
• networkx for manipulating graphs and networks

Major libraries you might use in data science include

• numba for just-in-time compilation of numpy-based code, which can be used to speed up computations and work around the lack of a compiler for Python.
• pandas for large dataset organization, similar to R’s data.frame and SQL tables
• statsmodels for statistics and data analysis, with API similar to R’s lm and glm
• seaborn for alternative plotting, with API similar to R’s ggplot
• PyTorch for deep learning.
• JAX for automatic differentiation and GPU acceleration, as well as deep learning

To use an external library that is not already installed with Python, you can install it using Python’s package manager, pip.

    pip install pandas

This will install whatever version of the package has been made available by the package maintainer on the Python Package Index. pip will also install all dependencies of the package, and update any existing installations or dependencies. After installing a package into Python, you usually must relaunch your Python kernel.

Dependency conflicts and virtual environments¶

When installing a single package like pandas, you will often find that pip installs a large number of other packages on which pandas depends, including numpy. However, if you later install another external package like scikit-learn, it may also depend on numpy. If scikit-learn requires the same version that pandas already installed, pip will not install a second copy, it will simply use the existing copy. Likewise, if the scikit-learn version is a newer release that requires a newer version of numpy, by the newer version is also compatible with pandas, then pip will go ahead and install the newer version, overwriting the older version.

An issue arises if the second package requires a different version of numpy than the first, and the two versions are not compatible. In that case, the package manager will run into conflicts. In large projects with many dependencies, managing these chains of dependencies can become quite difficult.

One solution to this problem is to use a virtual environment. A virtual environment is a self-contained Python installation, including a copy of the Python interpreter, a copy of the standard library, and a copy of all the external libraries you have installed. This allows you to have different versions of the same library installed for different projects, and it allows you to isolate the dependencies for one project from the rest of your projects. Most Python distributions include a package manager for creating and managing virtual environments, such as conda, mamba, uv, venv, and pyenv. Many of these environment managers also include their own alternative to pip for installing packages, which can manage conflicts.

Numpy¶

Numpy is technically an external library, but we will use it so often in scientific computing that it’s basically a built-in data type for our purposes. The key feature of numpy is the ndarray datatype, which represents tensors. Most functions for linear algebra and other scientific computing tasks are within numpy, and many share the same names and syntax as their MATLAB equivalents.

Arrays, indexing, broadcasting, vectorization¶

The logistic map¶

The logistic map is a recurrence relation that maps a single scalar value to another scalar value. It is often cited as an archetypal example of how complex, chaotic behaviour can arise from very simple nonlinear dynamical equations.

The map is defined by:

where the dynamical variable $x_n \in (0, 1)$ is a number between zero and one. The parameter $r$ represents a growth rate, and varies between zero and four.

For $r$ values between zero and one, the population always dies out. However, as $r$ increases past one, a bifurcation occurs at $r=3$, where the population begins to oscillate between two values. As $r$ increases further, period doubling occurs, leading to four values, then eight, and so on, until finally chaos occurs, and the population takes on seemingly random values between zero and one.

Runtime complexity and dynamical systems¶

A common question we might ask for an algorithm is: How does the runtime scale with the input size? What about the memory usage? These are sometimes called the asymptotic time and space complexities of an algorithm. The asymptotic complexity of an algorithm is independent of the hardware and software environment. It is a representationn-invariant mathematical property of the algorithm.

Notice that in the asymptotic regime the runtime scales linearly. Why is the runtime $\sim N$? We only “touch” each element of our solution array twice---once to calculate it, and once to update the next state. We don’t need to look back at all the previous states during each update, and so adding additional iterations doesn’t incur an unexpected cost. In our case, the memory usage is also $\sim N$ because we are storing the full, cumulative history in order to plot the trajectory. If we only wanted to plot the final state, we could get away with only storing the final state. In this case, the memory would be $\sim 1$.

Now that we have a fast implementation of the logistic map, we can explore its properties in more detail. We will create a bifurcation diagram, which shows the values of $x_n$ for each value of $r$. For many values of $r$, we simulate the population dynamics for a long time, discarding an initial transient. We then record all of the unique values of $x_n$ that we observe, and plot all of them as points versus that value of $r$. The bifurcation diagram is thus a one-to-many mapping from $r$ to multiple values of $x_n$.

Storing the unique values of $x_n$ is a little fuzzy since $x_n$ is a continuously-valued variable and so, in principle, every single $x_n$ is a different number to within machine precision, but if there are any signs of discretization, we can see them in the bifurcation diagram.

The logistic map follows the celebrated period-doubling route to chaos. As the parameter $r$ increases, the period of the oscillations doubles, until chaos occurs at $r \approx 3.57$.

Notice how the spacing in $r$ between the period-doubling bifurcations decreases as $r$ increases. The geometric rate of decrease in the spacing between bifurcations is a constant value called the Feigenbaum constant.

Slicing and indexing numpy arrays¶

Numpy arrays follow the same indexing rules as Python lists. We can slice them and index them using the same syntax. Below, we will use slicing to create a return map of the logistic map. This is a plot of $x_{n+1}$ vs. $x_n$.

We can also consider higher-order return maps, which plot $x_{n+k}$ vs. $x_n$. For example, for $k=2$

Strided slicing¶

We can use strided slicing to take every $k^{th}$ element of an array. To illustrate this effect, we wil first simulate the Logistic map for several values of the bifurcation parameter $r$. We will store all of these trajectories in a list.

We can see that these trajectories all look pretty different, but what happens if we resample them? We will use strided slicing to take every nth element of each array.

We can see a trivial self-similarity in the trajectories in the ordered phase: we can take every other element of a discrete period-4 trajectory, and it looks like a period-2 trajectory. Likewise, we can take every fourth element of a period-8 trajectory, and it looks like a period-2 trajectory. However, we can’t do this for the chaotic phase.

Recurrence plots¶

We can see that the logistic map exhibits markedly different dynamics in different regimes. The chaotic phase is erratic and non-repeating. However, if we encountered a time series from the chaotic phase, we might wonder whether it came from a chaotic deterministic system, or a stochastic system.

One way to address this question is to calculate the recurrence plot of the time series. The recurrence plot shows the distance between all points in the time series.

We have already calculated example trajectories for a logistic map in the ordered and chaotic regimes. We will also, for comparison, calculate a random time series.

Boolean Selection with numpy arrays¶

We can also index into the elements of Python arrays by using a Boolean array of the same shape. These will select out elements of the outer array wherever the Boolean array is True. Regardless of the shape of the outer array, the resulting array be a 1D array with length equal to the number of True elements in the Boolean array.

This function only differs from the original recurrence plot function by one line, but that line performs and elementwise comparison between the recurrence plot matrix and a float threshold. This returns a new matrix with the same shape as the original, but with type bool.

We can also use Boolean selection to select points from an array. For example, we can select all points from an array that are less than 0.5.

Finding the changepoints of a time series¶

Now, let’s use Boolean selection to select the points corresponding to when the logistic map crosses the point $x = 0.5$. We will first threshold the trajectory into a series of 0s and 1s corresponding to whether the trajectory is above or below 0.5. We then use np.diff, which calculates the difference between subsequent elements of the logistic map, in order to find places where the map crosses the threshold.

Interestingly, we observe an exponential waiting time distribution, in which the map occasionally takes a very long amount of time to cross from one side of the threshold to the other.

Fancy indexing¶

• Can pass lists of indices to modify

Broadcasting¶

We can slice numpy arrays just like lists, but in multiple dimensions. However, numpy arrays recognize a formal difference between singleton dimensions and empty dimensions.

In 1D, singleton dimensions distinguish between row and column vectors, and the same idea holds for higher-dimensional arrays. We will often exploit singleton dimensions in order to guide numpy regarding the linear algebra operations we want to perform.

We can lift an array to higher dimensions by adding singleton dimensions. This can be done by indexing into the array with np.newaxis or by using the None keyword.

Vectorization¶

Vectorization is the process of representing low-level array manipulations using high-level tensor operations in numpy. Usually, vectorization does not change the intrinsic runtime scaling of an algorithm, but it changes the prefactor. However, faster runtime can sometimes come at the expense of higher memory usage.

Let’s try finding the pairwise distance matrix between two sets of points in two-dimensions. We will do this in two ways: a naive way, and a vectorized way.

Mathematically, we have a vector $X$ of shape $(N, 2)$. We wish to compute the matrix $D$ of shape $(N, N)$, where $D_{ij} = ||X_i - X_j||_2$.

Now, let’s modify our approach to use a built-in numpy function, np.linalg.norm, which calculates the Euclidean norm of a vector. This function is internally vectorized, allowing it to process all dimensions a vector quickly.

Since we know the shape of our output matrix, we can also preallocate the memory for it using np.zeros. This is usually a good practice to avoid Python having to continually resize a list as we add elements to it. However, the speed-up is not as dramatic as in older languages like MATLAB, due to differences in the way that lists and arrays are implemented in Python.

Conceptually, we can think of vectorization as numpy processing all of the dimensions of the array at once, rather than looping over them one at a time.

We can do even better if we take the idea of vectorizing code a bit further. We think about the distinct “dimensions” present in our problem: the two dimensions of the output array (both have value 200), which each correspond to the first dimension of the input array, and the dimension corresponding to the second dimension of the input array, which has value 2. Symbolically, we start out with the input array $x_{id}$, where $i$ is the index over the rows ($N=200$), and $d$ is the index over the dimensions ($D=2$). We want to compute the output array $D_{ij}$, where $i$ is the index over the $N$ rows, and $j$ is the index over the $N$ columns.

We therefore expect that we need to perform an expansion first, to lift the dimensionality of our $(200, 2)$ input array to an intermediate with shape $(200, 1, 2)$. We can do this using numpy’s broadcasting rules. Since we want to find the difference between each element of the input array with itself, we need to do a broadcast subtraction. Symbolically, we can write this as:

The singleton dimension thus tells numpy that we want to split the first dimension of the input array into two separate dimensions for indexing purposes. In contrast, the last dimension of the input array is not split, because we want to keep the difference between the two points in the same dimension.

We next need to perform a contraction, which pulls the $(200, 200, 2)$ intermediate array back down to a $(200, 200)$ output array. Symbolically, we can write this as:

We can implement this operation using numpy’s built-in functions, which are optimized for linear algebra operations. Since we are computing the Euclidean norm, we can use np.linalg.norm to do this contraction in a single step. Other common contractions might include np.sum, np.mean

The scaling stays the same, but the prefactor changes¶

The runtime of both functions appears to increase faster than linearly with the number of points, $N$. This makes sense, because finding the distance between all pairs of $N$ objects requires $\sim \mathcal{O}(N^2)$ operations. But the prefactor (the scale of the y-axis) is much lower for the vectorized code. So vectorization offers a practical speedup, but not an intrinsic change in the complexity of the problem.Generally, we should think of runtime complexity as an intrinsic mathematical property of the algorithm, while the “prefactor” is a practical consideration that depends on the implementation and hardware.

Under the hood, in vectorized code, numpy still loops over the elements of the array, but defers this operation to C, a compiled language, rather than Python. So by expressing our program using higher-level mathematical objects, we gave numpy enough insight into the code to properly optimize it. Generally, in vectorization, we think of the compiled, behind-the-scenes code as being executed instantaneously or in parallel (even if it’s not actually in practice). Whenever you have an array operation that seems like it can be done in parallel, there’s likely a way to vectorize it.

Finding the fractal basins of attraction of the Mandelbrot set¶

The Mandelbrot dynamical system is defined by the following iteration:

where the complex number $z = x + i y$ encodes a two-dimensional spatial position. We can equivalently write this as a pair of real-valued update rules for the real and imaginary parts of $z$:

Depending on the value of $c$, the trajectory of $z_t$ may diverge to infinity, or converge to a finite value. The set of points $c$ for which the trajectory of $z_t$ does not diverge is known as the Mandelbrot set.

Implementing the Mandelbrot set¶

For our implementation of the Mandelbrot set, we could work directly with float arrays for $x$ and $y$. However, Python conveniently provides a built-in complex data type, which is a pair of floats representing the real and imaginary parts of a complex number. We will use this data type to represent the complex numbers in our Mandelbrot set.

Let’s start by plotting two trajectories of the Mandelbrot set that start from different initial conditions.

We are interested in finding the basins of attraction, which the subregions of the complex plane that converge to the same asymptotic value. We can do this by computing trajectories of the Mandelbrot set for a grid of initial conditions in the complex plane.

In our implementation, we iterate over a grid of initial conditions and then calculate the trajectory under the Mandelbrot map. If the absolute value of the trajectory exceeds a threshold, we assume that means that those initial conditions diverged. We therefore assign the initial condition to the basin of attraction corresponding to the threshold that was exceeded.

In our vectorized implementation, we will use the same algorithm, but we will iterate the entire initial condition array all at the same time. We will use Boolean indexing to assign different values, depending on whether a given array point is diverging or not.

Let’s do a timing experiment, where we see how the runtime of our two implementations changes as we increase the number of iterations.

Our numpy vectorized implementation is fast enough that we can zoom in and explore the fractal structure of the Mandelbrot set.