# Numpy: vectorisation and arrays

This short chapter introduces the essential ideas behind Numpy.
A longer and better version is the [Numpy fundamentals](https://numpy.org/devdocs/user/basics.html) in the Numpy User Guide, which you would do well to skim afterwards.
A proper (and very good) textbook for anyone already at home in Python – which is what we take ourselves to be at this point in the course – is [Nicolas Rougier's From Python to Numpy](https://www.labri.fr/perso/nrougier/from-python-to-numpy/); chapter 3 is the part that concerns us first.

In short: Numpy brings [array-oriented programming](https://en.wikipedia.org/wiki/Array_programming) (also called *vectorisation*) to Python.

Why is that worth having? Three reasons, and they hang together.

First, Numpy is above all an **interface to very fast libraries** for linear algebra and statistics – written in C, in Fortran and here and there in assembler, and optimised over decades. An `np.dot` calls code that nobody would want to rebuild in Python.

Second, this way of working allows **SIMD** (*single instruction, multiple data*): instead of looping over a list and handling each element on its own, you apply one operation to a whole vector. Modern processors can compute with several values in a single step – but only if the same thing is to be done to all of them. Branches (`if` tests inside the loop) ruin that. Hence this chapter's working rule: **write branch-free code, which is to say, vectorise.** Where a case distinction really is needed, you express it as a mask over the whole array rather than an `if` in the body of a loop; the section *Masks and fancy indexing* shows how.

Third, this way of thinking carries over. PyTorch, the standard library for neural networks, works in almost the same way – there the object is called a `tensor` rather than an `ndarray`, and it lives either in main memory or on the graphics card. Whoever has learned to vectorise here can carry on doing it there; the course *Introduction to Deep Learning* assumes precisely that.

We use Numpy version 2 or higher, for which Python 3.11 or above is recommended. Most of the code in these notes should run under Numpy 1.26 as well, though. [NEP 29](https://numpy.org/neps/nep-0029-deprecation_policy.html) has more on version compatibility.

## Arrays
<!--
On probability with Numpy, this one is a good fit (in German):
https://www.python-kurs.eu/python_numpy_wahrscheinlichkeit.php
-->
Numeric Python (Numpy) is usually imported under the abbreviation `np`:

```python
import numpy as np
np.__version__
```

```python
sample = np.random.random(10)
print(sample)
print(type(sample))
sample
```

Numpy contains a great many useful helper functions, which in turn operate on Numpy's own data structure, the array.

```python
sample = np.random.randint(low=0, high=10, size=5000)
print("mean    =", np.mean(sample),   "\nexpected=", 9/2)
# 99/12 = (n²-1)/12 with n=10, the variance of the discrete uniform distribution
# (derived in the chapter "The most important discrete distributions"):
print("std²    =", np.std(sample)**2, "\nvariance=", 99/12)
```

There is an `arange` function too, but in Numpy it does not produce range objects:

```python
myRange = np.arange(1, 7)
print(myRange, type(myRange), myRange.dtype)
myRange
```

Unlike a Python list, the Numpy array `ndarray` carries one fixed data type that all its elements share. That can be any of the Numpy data types (dtypes): `double`, for instance (compatible with the Python `float`), `float32` (on most platforms compatible with the C `float`) or `long` (compatible with the Python `int`).

```python
import timeit
ordinary_list = [1,2,3,4,5,6,5,4,3,2,1]*10
def sort_array(dtype):
    a = np.array(ordinary_list, dtype)
    a.sort()
print(timeit.timeit(lambda : sort_array(np.byte)))
print(timeit.timeit(lambda : sort_array(np.float64)))
```

Because the data type is known exactly, Numpy can use algorithms optimised for it and implemented directly in C.

The documentation lists the other `dtype`s and how to use them:
* ["Structured Arrays"](https://numpy.org/doc/stable/user/basics.rec.html)
* ["Data types" in the "Numpy fundamentals"](https://numpy.org/doc/stable/user/basics.types.html)
* ["Scalars"](https://numpy.org/doc/stable/reference/arrays.scalars.html)
* ["Data type objects (dtype)"](https://numpy.org/doc/stable/reference/arrays.dtypes.html)
* ["numpy.dtype"](https://numpy.org/doc/stable/reference/generated/numpy.dtype.html)

Arrays can be created in more ways than `arange` or converting a Python sequence: `ones`, for instance, gives you an array filled with $1$, and `zeros` one filled with $0$.

The array is called `ndarray` because that stands for "$n$-dimensional array".
If you wanted to store a matrix in plain Python, you would do it as a list of the row vectors (or of the column vectors), something like

```python
matrix = [[1,0,0], [0,1,0], [0,0,1]]  # identity matrix
print(matrix, type(matrix))
shamtrix = matrix[:]   # a "copy" by slicing
print(shamtrix[0][0])
shamtrix[0][0] = 0     # we change the top left entry
print(matrix[0][0])
```

This exposes a problem: Python treats our matrix as a list – which is, after all, how we wrote it down – and so copying the list does not copy its contents (the row vectors), only the pointers to them.

```python
npmatrix = np.identity(3, int)  # identity matrix
print(npmatrix, type(npmatrix))
npshamtrix = npmatrix[:]
print(npshamtrix[0][0])
npshamtrix[0][0] = 2     # we change the top left entry
print(npmatrix[0][0])
# With a real copy this would not have happened:
real_copy = npmatrix.copy()
real_copy[0][0] = 1
assert npmatrix[0][0] != 1
```

Note that the Numpy array behaves exactly like our Python list of lists: we index and slice it the same way, and copying it via the slicing syntax has the same problem.

The `shape` parameter tells us what shape our Numpy array has. It is a $d$-tuple, where $d$ is the dimension. A matrix is $2$-dimensional, a vector $1$-dimensional and a scalar $0$-dimensional.

```python
print(npmatrix.shape)
print(npmatrix[0], npmatrix[0].shape)
print(npmatrix[0][0], npmatrix[0][0].shape)
```

Numpy lets you slice rather more finely.
The general syntax is `[start:stop:step, ..]`, where the comma moves you from one axis to the next. A two-dimensional array has two axes: axis $0$ runs from top to bottom and axis $1$ from left to right. The "step" works with Python lists too; indexing along several axes is a Numpy speciality.

```python
numbers = list(range(1,10))
npmatrix = np.array(numbers)
print("not in shape:", npmatrix)
npmatrix.shape = (3,3)
print("in much better shape:\n"+ str(npmatrix))
print("row vector, row 0:", npmatrix[0])
print("column vector, column 0:", npmatrix[:,0])
print("columns 1-2:\n"+ str(npmatrix[:,1:]))
print("all rows, step 2\n"+ str(npmatrix[0::2]))
```

In classic Python we would reach an entry of a matrix with `matrix[x][y]`, and that works in Numpy too. But it first creates a further array object `matrix[x]` (a view, see below) and only then calls `[y]` on it. In principle, then, Numpy's own `[x,y]` is the more efficient choice.

When slicing, Numpy deliberately makes no copies, only a different view of the same region of memory (which is also why `[:]` behaved as it did above). `np.may_share_memory` checks whether two arrays point into the same region of memory. A positive answer by no means says that the arrays depend on one another – the first and the second column of a matrix also point into the same region of memory, namely the whole matrix. Change one of the two vectors and the other stays as it was – but the whole matrix changes along with it.

```python
npmatrix = np.array(list(range(1,10))).reshape(3,3)
candidates = (npmatrix[0::2], npmatrix[1])
print("May share memory (but actually don't):",
      candidates, np.may_share_memory(*candidates))
print(type([0,0,0]))     # before the assignment
npmatrix[1] = [0,0,0]
print(type(npmatrix[1])) # after the assignment
print(npmatrix)          # the whole matrix has changed with it
```

For the details of slicing and indexing, see [the documentation](https://numpy.org/doc/stable/user/basics.indexing.html).

## Broadcasting

For Python lists, the addition operator means concatenation, which also gives multiplying a list by a scalar a well-defined meaning, while multiplying two lists is undefined.
Numpy arrays support a good many more arithmetic operations:

```python
E = np.identity(3, int)
print(E, "= E")
A = np.ones((3,3), int)
A[0] = [0,0,0]
print(A, "= A")
print(E + A, "= E + A")
print(E * A, "= EA")
print((E+A)**2, "= (E+A)(E+A)")
```

Look closely and you will notice that these are not matrix multiplications:
the calculations have simply been carried out element by element, because that is how the operations on arrays are defined.
This one is particularly treacherous:

```python
v = np.ones(3, int)
print(v, "= v")
print(A*v, "= A*v (but not the matrix product)")
```

For explicit matrix algebra, people used to reach for the Numpy data type `matrix`, but it is marked as `deprecated` and will disappear in future Numpy versions. These days you use `np.matmul` or the infix operator `@`.

```python
print(np.matmul(A,v))
print(A@v)
```

For many problems it is a great help to leave matrix algebra aside and carry out the arithmetic element by element.
*Broadcasting* is the mechanism that makes such element-wise arithmetic rather more convenient.
The operation ($n \times n$ array) * ($n$-vector), for instance, is interpreted automatically by copying the $n$-vector $n$ times, so that each row of the left-hand array is multiplied (element by element) by the vector.

It really is worth skimming [the documentation](https://numpy.org/doc/stable/user/basics.broadcasting.html) on this once.

A *tensor* in PyTorch, the library of the course *Introduction to Deep Learning*, is essentially a Numpy array with two additions: a device assignment (CPU or GPU) and automatic differentiation. `shape`, axes (`axis`) and broadcasting work there exactly as they do in Numpy, which is why that course assumes this chapter.

## Masks and fancy indexing

A comparison such as `a > 2` is evaluated element by element and yields an array of `True`/`False`: a **mask**. A mask can serve as an index: you get exactly the entries where it says `True`, and you can assign to those entries as well. This is the tool we filter data with in the next chapter ("all passengers under 13", "all rows with no missing value"). `np.where(condition, a, b)` chooses element by element between two arrays: where the condition holds, the entry from `a`, otherwise the one from `b`. That is how a per-element `if` test becomes a single array operation:

```python
a = np.array([3, 1, 4, 1, 5, 9])
mask = a > 2
print(mask)                        # mask: array of True/False
print(a[mask])                     # only the entries where it says True
print(a[(a > 2) & (a < 9)])        # combine conditions with & (and), | (or), ~ (not),
                                   # not with and/or; every condition in brackets!
a[a == 1] = 0                      # assigning to a mask
print(a)
print(a[[0, 2, 4]])                # "fancy indexing": a list of positions as the index
print(np.isnan(np.array([1.0, np.nan, 3.0])))   # masks for missing values
print(np.where(a > 2, a, 0))       # np.where: where a > 2 the value from a, otherwise 0 (no if needed)
```

## Probability vectors

In the chapter on random variables, a distribution on $n$ values was a vector $p \in \Delta_n$: non-negative entries that sum to $1$. In Numpy that is simply an array divided by its own sum, and a sample from this distribution is a call to `choice` with the argument `p`:

```python
p = np.array([1, 1, 1, 1, 1, 3], dtype=float)   # loaded die: the 6 is three times as likely
p = p / p.sum()                                   # normalise: p now lies in the simplex Δ_6
print(p, p.sum())

rng = np.random.default_rng(1)
rolls = rng.choice(np.arange(1, 7), size=10000, p=p)   # sample from the pmf p
print(np.bincount(rolls)[1:] / len(rolls))             # relative frequencies, roughly p
```

## ufuncs

`ufunc` stands for "universal function" and denotes a function that can run vectorised on Numpy arrays.

Every `ufunc` supports a `reduce`, for instance, which applies the `ufunc` along one axis of the array to the smaller arrays that result.

To write a `ufunc` of your own you have to [program in C](https://numpy.org/doc/stable/user/c-info.ufunc-tutorial.html) or else [wrap a Python function](https://numpy.org/doc/stable/reference/generated/numpy.frompyfunc.html).

It is worth taking a quick look at all the `ufunc`s that are already defined:
[Available ufuncs](https://numpy.org/doc/stable/reference/ufuncs.html#available-ufuncs)

## What do you need to know?

Numpy and the other libraries of scientific computing in Python are large; nobody can know every function.
For these notes, a pass over the code cells of all the notebooks (`grep -ohE "np\.[a-z_]+(\.[a-z_]+)?" *.ipynb | sort -u`, tidied up by hand) yields the complete list below, as of September 2026. You do not have to learn it by heart, but it does give a clear upper bound: more than this you definitely do not need to know! The functions in **bold** are the ones you really should have in your head, because you type them constantly; the rest is a reference work. Bear in mind, too, that the list covers the whole of these notes, and that much of it we have not met yet.

For day-to-day work a cheat sheet helps more than a list: the [Numpy reference](https://numpy.org/doc/stable/reference/index.html) has an overview sorted by topic, and for getting started there is the [compact beginners' chapter of the Numpy documentation](https://numpy.org/doc/stable/user/absolute_beginners.html). Print it out, keep it beside you, use it – learning it by heart gets you nowhere here.


 * **abs**
 * **all**
 * **allclose**
 * append
 * **arange**
 * **argmin**
 * **array**
 * atleast_1d
 * bincount
 * byte
 * c_
 * column_stack
 * **concatenate**
 * **copy**
 * **corrcoef**
 * cos
 * **cov**
 * **cumsum**
 * diag
 * **diff**
 * **dot**
 * dstack
 * equal
 * **exp**
 * **eye**
 * float64
 * **full**
 * genfromtxt
 * heaviside
 * **hstack**
 * identity
 * **isclose**
 * **isnan**
 * linalg.det
 * linalg.eig
 * linalg.eigh
 * **linalg.inv**
 * **linalg.norm**
 * linalg.svd
 * **linspace**
 * **log**
 * log2
 * logspace
 * **matmul**
 * **max**
 * **maximum**
 * may_share_memory
 * **mean**
 * **median**
 * **meshgrid**
 * mgrid
 * **min**
 * **minimum**
 * **nan**
 * nansum
 * **ones**
 * ones_like
 * outer
 * **pi**
 * power
 * ptp
 * **random.choice**
 * **random.default_rng**
 * random.multivariate_normal
 * random.permutation
 * random.rand
 * random.randint
 * random.randn
 * random.random
 * **random.seed**
 * random.uniform
 * **round**
 * searchsorted
 * set_printoptions
 * sin
 * size
 * **sort**
 * **sqrt**
 * **square**
 * **stack**
 * **std**
 * **sum**
 * **unique**
 * **var**
 * **vstack**
 * **where**
 * **zeros**
 * zeros_like