NumPy 1.20 has been released

NumPy is a Python library that adds an array data type to the language, along with providing operators appropriate to working on arrays and matrices. By wrapping fast Fortran and C numerical routines, NumPy allows Python programmers to write performant code in what is normally a relatively slow language. NumPy 1.20.0 was announced on January 30, in what its developers describe as the largest release in the history of the project. That makes for a good opportunity to show a little bit about what NumPy is, how to use it, and to describe what's new in the release.

What is NumPy?

NumPy has been a big part of what makes scientific Python possible. It is the library at the core of SciPy and most other science or engineering ecosystems for Python. It allows researchers to use a scripting language with a friendly syntax for the analysis of data from simulations and experiments. It can even allow them to run many types of simulations in Python itself. NumPy is largely responsible for popularizing the idea that Python can be used for the kind of work that, previously, had been the almost exclusive domain of Fortran. This has had a cascading effect, leading to the wider appreciation of open-source tools and libraries by the scientific community, and to the emergence of new modes of sharing results; in particular, Jupyter notebooks.

NumPy is released under the BSD 3-Clause License and is developed on GitHub.

NumPy adds a new data type to Python: the multidimensional ndarray. This a container, like a Python list, but with some crucial differences. A NumPy array is usually homogeneous; while the elements of a list can be of various types, an ndarray will, typically, only contain a single, simple type, such as integers, strings, or floats. However, these arrays can instead contain arbitrary Python objects (i.e. descendants of object). This means that the elements will, for simple data types, all occupy the same amount of space in memory. The elements of an ndarray are laid out contiguously in memory, whereas there is no such guarantee for a list. In this way, they are similar to Fortran arrays. These properties of NumPy arrays are essential for efficiency because the location of each element can be directly calculated.

Beyond just adding efficient arrays, NumPy also overloads arithmetic operators to act element-wise on the arrays. This allows the Python programmer to express computations concisely, operating on arrays as units, in many cases avoiding the need to use loops. This does not turn Python into a full-blown array language such as APL, but adds to it a syntax similar to that incorporated into Fortran 90 for array operations.

Before delving into more detailed examples in the next section, let's take a brief look at how this works. Suppose we have a Python program with two lists that we want to add together, element by element, to make a new list. Simple code to do so might look something like this:

    >>> a = [1, 2, 3]
    >>> b = [4, 5, 6]
    >>> c = []

    >>> for i in range(len(a)):
    ... c.append(a[i] + b[i])

After this code is run, c will be [5, 7, 9].

The equivalent code using NumPy will look like this:

    >>> import numpy as np     # the conventional NumPy import
    
    >>> a = np.array([1, 2, 3])
    >>> b = np.array([4, 5, 6])

    >>> c = a + b

The library overloads the "+" operator to act element-wise on the NumPy array type, so no looping is required. The result, c, will not be a list but another ndarray:

    >>> c
    array([5, 7, 9])

    >>> type(c)
    numpy.ndarray

The same syntax can be used for arrays of any dimension. Element-wise multiplication of two matrices is accomplished in the obvious way:

    >>> a = np.array([[1, 2], [3, 4]])
    >>> b = np.array([[2, 2], [2, 2]])  

    >>> a * b

    array([[2, 4],
           [6, 8]])

But what if we want matrix multiplication, rather than element-wise multiplication? NumPy also implements the binary operator "@" for that:

    >>> a @ b

    array([[ 6,  6],
           [14, 14]])

    >>> b @ a

    array([[ 8, 12],
           [ 8, 12]])

The above examples show how NumPy extends Python arithmetic to operate on entire arrays, but the library offers far more than this. There is a complete set of facilities for creating and manipulating arrays, including transposing, rotating, inverting, joining, and selecting sub-arrays; a sub-library for Fourier transforms; a set of linear algebra routines including such things as computing eigenvalues, determinants, and solving systems of linear equations; and more. NumPy also includes a subsystem called F2PY for creating interfaces between Python scripts and Fortran subroutines.

Much of this capability, especially in the realm of linear algebra, is provided by the widely-used BLAS and LAPACK libraries, for which NumPy acts as a convenient wrapper for the Python programmer. The user has the choice of using the versions supplied with NumPy or installing [147-page PDF] a hardware-specialized version for maximum efficiency.

NumPy performs admirably when used with standard array calculations that fit the patterns that it is designed to optimize, such as array multiplication or element-wise arithmetic. Its advantages are limited to such calculations, however, and this can be considered as the boundary that delimits the practical usefulness of Python for scientific computing.

More complex algorithms that cannot be expressed as a sequence of simple array operations, and therefore require inner loops and conditionals, will run at the speed of the Python interpreter, and not scale to large problems. The user then has several choices: to turn to a language designed for high performance computing, such as Fortran or Julia, to resort to alternative implementations of Python, such as Numba, PyPy, or Cython, or to develop in two languages, using Python as glue and a fast language for the numerically intensive parts. Here is a comparison of some of these approaches, which also serves as a brief tutorial, with examples.

Examples

In order to get the recent release of NumPy, the easiest route is probably to use pip to get it from the Python Package Index (PyPI). As usual with Python, the standard recommendation is to do this installation, and the installation of anything else to be used with numpy, such as IPython, in a virtual environment.

In my recent article about SciPy, I demonstrated the use of its image-processing routines to show some of what it can do. Under the hood, all of these image transformations are accomplished by converting the image to a matrix of numbers and using NumPy to transform the matrix in various ways. It will be instructive to carry out some similar image transformations now, not by using the image-processing library, but by using NumPy array manipulations directly. This will demonstrate some of the things that can be done with NumPy and its arrays in a way that makes the operations more intuitive because of the direct visual results.

First we need to import numpy and the plotting library (needed to display the results as images) as well as read in an image file and store it as a NumPy array. All the examples in this article use NumPy 1.20.0 and Python 3.9.1.

    >>> import numpy as np
    >>> import matplotlib.pyplot as plt

    >>> img = plt.imread("bacall.jpg")

Now img will be a NumPy array. We can check this, and also get its dimensions:

    >>> type(img)
    numpy.ndarray

    >>> np.shape(img)
    (550, 440)

The image has dimensions 550 by 440. We can prepare this image using the following code:

    >>> plt.imshow(img, cmap='gray')

That tells plt that we would like the intensities mapped to a grayscale palette. Now we must either call plt.savefig() to save the image in a file or plt.show() to display it on the screen. We'll omit these calls, from here on, for brevity, just showing the matrix transformations.

Looking at the effect on images is a good way to check whether our mental model of certain matrix operations is correct. The image below shows a collection of operations done on the original image. NumPy has a built-in transpose() operator, that switches axes. "Transposed" is the result of np.transpose(img). Certain operations are easier to show than to explain in words. See the "Rolled" image for the result of np.roll(img, 200). Combining flip(), which does what it sounds like, with hstack(), which concatenates arrays horizontally, with np.hstack((img, np.flip(img))) generates the "Flipped" image; note that hstack() takes a tuple argument, thus the "extra" parentheses.

We could implement more sophisticated image processing, such a blurring, filtering, edge detection, and so on, by combinations of Fourier transforms, numerical derivatives, and other straightforward operations on the image matrix using NumPy.

The new release

There are several aspects to the new release of NumPy that will be significant for current users; some of the changes also give indications about NumPy's intentions for the future. First of all, NumPy 1.20 requires at least Python 3.7; no earlier version of Python will work. Beyond that, removing compatibility with Python 2 allowed additional code cleanup and elimination of what the release notes called "technical debt".

One class of techniques that NumPy uses to speed up Python code is to perform certain calculations in parallel, when the opportunity is afforded by the hardware that the code is running on. There are many forms of parallel processing. The general category that NumPy implements is single instruction, multiple data (SIMD) parallelism. There are two kinds of SIMD calculation, and both involve performing arithmetic on different elements of an array simultaneously. The first kind is vectorization, where a single CPU (or core) carries out an arithmetic operation on a small number of elements, typically four or eight, in parallel, continuing until the entire array is processed. The second kind of SIMD parallelism divides the array among a set of processors, from the small handful found on a laptop to the thousands of processors on a supercomputer. This is the form of parallelism provided by GPU array processors.

NumPy takes advantage of both types of SIMD parallelism when possible. Multicore, distributed parallel computation is provided by the BLAS routines wrapped by NumPy; whenever a NumPy operation calls a routine from this library, it takes over the parallelization of the calculation, and will utilize all of the processing units available, transparently to the user.

For example, matrix multiplication using the @ operator calls, on my machine, an OpenBLAS routine that uses all the cores on my laptop. If I had installed NumPy using Anaconda, rather than pip, the same Python code would call out to the Intel MKL library, which also implements distributed SIMD parallelism.

The typical user need not be aware of any of this, but sophisticated users have always been able to choose the best fit for their target hardware from among several available linear algebra libraries, and compile a custom version of NumPy that links to it. The new release of NumPy extends similar flexibility and control into the realm of vectorization. The main advance is, again, something that most people need not be aware of, and the only visible change to most users will be a possible increase in performance. This increased performance is due to the new ability of NumPy to detect, at runtime, which vectorization instructions are provided by the CPU, and choose a code path that takes advantage of them. The same line of Python will now execute different code, depending on what NumPy can detect about the hardware. One side effect of this is that NumPy is a bit bigger, because the compiled library now contains code paths for different CPU architectures.

As part of this initiative, there are several new build arguments for use in compiling custom versions of NumPy. These allow the user to select from among various vectorization optimizations, to help in tuning NumPy to a particular architecture, or even to a particular computation pattern. To aid this even further, there is now a #define called NPY_DISABLE_OPTIMIZATION that can be used to turn off optimizations for particular sections of the NumPy source code.

Another significant new feature in NumPy 1.20 is the type annotations. Optional static type annotations were added to Python with PEP 484, an addition to the language that was not without controversy. These type hints can be used by linters, type checkers, and development environments, especially to assist with code completion and refactoring.

NumPy now embraces PEP 484 type annotations in two ways. The NumPy library code is now annotated with type information. This had been discussed before, but dropping Python 2 support in the latest release of NumPy cleared the way for adding type annotations without the awkwardness of supporting Python versions where they can't be used. An examination of the files in the NumPy package that I installed shows that annotations have been implemented as stub files.

The second new feature related to type annotations is the addition of new types, defined and implemented in a new NumPy submodule called numpy.typing. The two new types mentioned in the documentation are ArrayLike and DTypeLike. The rationale is to enable the use of type checkers to help the user avoid patterns that, while legal in NumPy, would be inefficient. While using arbitrary objects in an ndarray is legal, it does not generally work with BLAS and other libraries that expect simple values. The new static types can be used to warn the user if they create such data structures.

In addition to the two new types mentioned in the documentation, the typing module contains definitions for a handful of others, such as FloatLike, ComplexLike, and ShapeLike. NumPy type annotations are still at an early phase and users can expect them to evolve significantly.

Another addition is a new keyword argument, where, in the numpy.all() and numpy.any() functions. These functions return True or False if any or all elements of an array satisfy a condition. For example, suppose we define an array containing some random integers this way:

    >>> ri = np.random.randint(0, 20, 8)
    >>> ri
    array([ 9,  7,  6,  2, 14, 13, 19, 16])

The random.randint() function takes a minimum, maximum, and an output shape, and returns an array of uniformly distributed random integers in the given range. We can ask if any or all of the integers are less than three:

    >>> np.any(ri < 3)
    True
    >>> np.all(ri < 3)
    False

Using the new argument, we can limit which elements are looked at in the test. The where argument must be set to an array of Boolean values, which is used as a mask to pick out the elements to be tested. NumPy makes it easy to create such logical masks with a simple Boolean statement. Here we create a mask that maps which elements are even, using the ndarray modulus operator:

    >>> ri%2 == 0
    array([False, False, True, True, True, False, False, True])

Now we can ask "are any of the even (or odd) numbers in the array less than three?":

    >>> np.any(ri < 3, where=ri%2 == 0)
    True
    >>> np.any(ri < 3, where=ri%2 == 1)
    False

This new argument leads to more succinct code, and more opportunities for parallelization in certain circumstances. It can also be applied to three statistical functions: mean(), std(), and var(), for the mean, standard deviation, and variance, as, for example, a concise way to eliminate outliers from the computation of these distribution parameters.

The new function sliding_window_view() provides a convenient way to generate a moving window over an array for doing such things as creating moving averages. This is quite convenient, since the alternatives using index arithmetic can lead to messy code. A simple example will serve better to show how this works than an attempt to explain it in words:

    >>> a = np.array([1, 2, 3, 4, 5, 6, 7])
    >>> np.lib.stride_tricks.sliding_window_view(a, 3)
    array([[1, 2, 3],
	   [2, 3, 4],
	   [3, 4, 5],
	   [4, 5, 6],
	   [5, 6, 7]])

As the example shows, this lives in the stride_tricks module inside the NumPy lib module. The documentation cautions that this function can be quite slow; it is intended for prototyping or for small arrays.

Many of the changes in the new release are deprecations involving moving functions, either from the top-level numpy namespace into submodules, or up from submodules into the top level, and the spelling of data types. Examples are the financial functions (present and future value, rates of return, etc.), now moved into their own numpy-financial package; the functions in numpy.dual have now moved into the top level.

The spellings of some data types have also been changed: type codes such as Uint32, that begin with upper-case letters, are now spelled uint32; the former style is deprecated. Certain data types that were denoted as, for example, numpy.int are now to be called simply int; see this table for a complete list. There were some possibly confusing notations for complex number data types, previously deprecated, whose deprecation periods have expired with the new release. These were data types such as Complex64, where the "64" referred to the size of the real or imaginary part separately. Users must now use complex128 instead of Complex64.

Historically, deprecations in NumPy have had long sentences in purgatory, at least one that lasted 14 years. Nevertheless, the reorganization that is part of this large release might be a good opportunity for cleaning up deprecated references in existing code bases.

The scale of the work devoted to the new release of NumPy, especially in the area of SIMD optimization, shows that its large community of developers see a continuing future for Python as a partner in scientific research and engineering calculation. There are many excellent alternatives for high-performance scientific numerics, from Julia to Fortran. Part of what Python brings to the table is the network effect of its large user base, and the existence of many solid scientific libraries. The consequence is that there is a good chance that a computational scientist today will wind up programming using NumPy. It's clear that this infrastructure is on solid ground, has a lot of talent behind it, and will only continue to improve.