2.2 Introduction to NumPy


: 40 minutes

We will use NumPy throughout this course to manipulate tensors—called ndarrays (n-dimensional arrays, or simply arrays) in NumPy. NumPy, short for Numerical Python, is the foundational package for scientific computing with Python. Pandas, scikit-learn, SciPy, and matplotlib are all built on top of it, so time spent here pays interest for the rest of the semester.

NumPy arrays provide fast, vectorized arithmetic and flexible broadcasting. Vectorized operations typically run one to two orders of magnitude faster than the equivalent Python for loop, because the loop happens in compiled C over a contiguous block of memory rather than in the interpreter over a list of boxed objects.

Import NumPy into your workspace in the following way:

NumPy Arrays

NumPy offers its N-dimensional array object, or ndarray, to represent a tensor and to facilitate efficient operations with them. Throughout this course, an array will refer to a NumPy N-dimensional array object, and a tensor its mathematical counterpart.

Recall that the matrix \bold{A}=\begin{bmatrix}1.5 & -0.1 & 3\\ 0 & -3 & 6.5\end{bmatrix} is a 2nd-order tensor with shape (2, 3).

To give you a taste of how NumPy works, we create an array named data to represent \bold{A}:

One can then perform mathematical operations with data. Note that no loop appears anywhere:

dtype

An array is a container for homogeneous data. Unlike a Python list, every element of an array has the same type, and that type is fixed when the array is created. NumPy offers many data types. For most purposes in this course we will meet float64 (double-precision floating point), int64 (long integers), and bool.

This homogeneity is not a limitation to be worked around; it is the reason arrays are fast. Because every element occupies the same number of bytes, NumPy stores the whole array as one contiguous block and can hand it to compiled code without inspecting anything element by element.

It also means dtype is sticky, and that has consequences:

Assigning 2.7 into an integer array does not promote the array to floating point; it truncates the value to 2 and moves on, without a word of complaint. If you need decimals, say so at creation time.

shape

Every array has a shape: the shape of the tensor it represents. It is a tuple giving the size along each axis.

These four attributes are the first thing to print when something is not behaving. Between them they answer nearly every “why did that happen?” question in this module.

reshape

The reshape method returns a new view of the same data laid out under a different shape.

Two facts about reshape are worth committing to memory.

First, it must preserve size. You cannot reshape 6 elements into a (2, 4) array, and NumPy will say so:

Second, you may leave one axis as -1 and let NumPy solve for it from the size. This is by far the most common way reshape is written in practice, because it stops you hard-coding a number that might change:

Two -1s are ambiguous and are rejected. And note the order in which elements are refilled: NumPy walks the last axis fastest (so-called row-major or C order), which is why np.arange(6).reshape(2, 3) puts 0, 1, 2 in the first row rather than the first column.

Array Creation

One can use the function np.array to very conveniently create a new array from an iterable Python object like a list, list of lists, and dictionary.

From a list

As an optional argument, you can pass the desired data type.

From a list of lists

WarningAutomatic data type

You will have noticed, on running the chunk above, that NumPy printed a decimal point after every value—even the whole numbers. Because the input list contains a floating-point value (1.2), NumPy picks a dtype wide enough to hold everything, namely float64, and converts the integers to match. Confirm it yourself by running B.dtype.

The general rule: NumPy infers the narrowest dtype that loses no information. It is usually best to let it do so rather than to force a dtype by hand.

Forcing a narrow dtype is a genuine foot-gun, because the values have to fit. An int8 holds only the range -128 to 127:

The value 152 does not fit in an int8, so it wraps around to -104. No error, no warning, just a wrong number sitting quietly in your data. Narrow integer dtypes save memory and cost correctness; unless you have a specific reason, take the default.

numpy.arange

Very much like Python’s built-in range(), numpy.arange creates evenly spaced values. The general form is np.arange(start, stop, step), and the interval is half-open: start is included, stop is not.

Unlike range(), the arguments need not be integers:

Note that start and step are optional, defaulting to 0 and 1.

Cautionarange with a fractional step is unreliable

How many elements does np.arange(0, 1, 0.1) have? The obvious answer is 10, and here it happens to be right—but the count is computed as \lceil(\textrm{stop}-\textrm{start})/\textrm{step}\rceil in floating-point arithmetic, and 0.1 is not exactly representable in binary. With a slightly different step the rounding can tip the other way and you get one element more, or fewer, than you expected.

When the step is fractional and the count matters, use linspace instead. Reserve arange for integer steps.

numpy.linspace

numpy.linspace(start, stop, num) asks the opposite question: instead of specifying the spacing and letting the count fall out, you specify the count and let the spacing fall out. Both endpoints are included by default.

Read that as “6 points from 1 to 4 inclusive”, so the gap is (4-1)/(6-1) = 0.6. The denominator is num - 1, not num, because you are counting intervals between points—an off-by-one that catches everyone once.

numpy.eye

This function creates a 2D identity array (or matrix) of a given dimension.

numpy.diag

numpy.diag can be used to create a 2D array (or matrix) with a prescribed vector as its main diagonal as the shown in the following:

A second argument k shifts the diagonal away from the main one; positive k moves it up and to the right, negative k down and to the left. Note that the result grows to accommodate the offset—four values on the first superdiagonal need a 5\times 5 matrix.

numpy.diag is its own inverse, in the following sense: given a 2D array it goes the other way and extracts the diagonal as a vector.

So np.diag decides what to do from the order of its input, which is convenient in a script and confusing in a traceback. If in doubt, print the shape.

numpy.zeros and numpy.ones

For a desired shape an array can be created filled with either zeros or ones. Note that the shape is passed as a single tuple, not as separate arguments:

Random Number Generator

NumPy provides pseudo-random number generation through its numpy.random module, which you should prefer to Python’s built-in random because it generates whole arrays at once.

The modern interface is a generator object, built by default_rng. You create one, seed it, and then call methods on it:

Common methods of a numpy.random.Generator
Method Description
random Draw samples uniformly from [0, 1)
uniform Draw samples uniformly from a given [\textrm{low}, \textrm{high}) range
integers Draw random integers from a given low-to-high range
standard_normal Draw from a normal distribution with mean 0 and standard deviation 1
normal Draw from a normal (Gaussian) distribution with a given mean and standard deviation
binomial Draw from a binomial distribution
beta Draw from a beta distribution
chisquare Draw from a chi-square distribution
gamma Draw from a gamma distribution
permutation Return a randomly permuted copy of a sequence
shuffle Randomly permute a sequence in place
choice Draw a random sample from a given array, with or without replacement
ImportantSeed once, at the top

default_rng(42) fixes the stream, which is what makes an analysis reproducible—and what makes these lecture notes print the same numbers for you as for me.

The idiom is to create one generator and reuse it:

Re-seeding before every call is a classic error, because it resets the stream and hands you the same “random” numbers over and over:

For more, read the documentation.

Exercises

Note Counting without running

How many elements does np.arange(2, 11, 0.3) contain? Work it out on paper from the half-open interval [2, 11) and the step 0.3, then assign the integer to answer.

Do not compute it with len(np.arange(...))—the point is to be able to predict it.

answer = 30

The span is 11-2 = 9 and the step is 0.3, so 9/0.3 = 30 steps fit, giving 30 values from 2.0 up to 10.7. The endpoint 11 is excluded.

Note Solving for the missing axis

flat holds 36 values. Reshape it into a 2-D array with exactly 9 columns, without typing the number of rows yourself—make NumPy work it out.

answer = flat.reshape(-1, 9)

flat.reshape(4, 9) gives the same array, but the -1 form states the intent (“nine columns”) and survives a change in the number of observations.

Note arange versus linspace

You need an array of exactly 5 values, evenly spaced, beginning at 0 and ending at 1, with both endpoints included. Which call does the job?

np.linspace(0, 1, 5) gives [0., 0.25, 0.5, 0.75, 1.]—five values, both endpoints present.

The distinction between the two functions is worth internalising:

  • arange(start, stop, step) fixes the spacing and lets the count fall out, over a half-open interval. np.arange(0, 1, 0.2) therefore stops at 0.8 and never reaches 1.
  • linspace(start, stop, num) fixes the count and lets the spacing fall out, and includes stop by default. The spacing is (\textrm{stop}-\textrm{start})/(\textrm{num}-1), since 5 points bound 4 intervals—hence 0.25 here, and hence linspace(0, 1, 4) giving a spacing of 1/3 rather than the requested five values.

There is a second reason to prefer linspace whenever the step is fractional: arange computes its length in floating-point arithmetic, and a step such as 0.1 is not exactly representable in binary, so the count can come out one larger or one smaller than you predicted.

Note Homogeneous Data

Consider the following, which runs without any error or warning.

counts = np.array([10, 20, 30])
counts[1] = 25.9
print(counts)

What is printed, and why?

An ndarray is homogeneous: its dtype is fixed when the array is created and assignment never changes it. Since [10, 20, 30] contains only integers, the inferred dtype is an integer type, and storing 25.9 into it casts toward zero, giving 25—not 26, and not a widened float array.

This is precisely the sort of bug this module is about: no exception, no warning, just a wrong number. A Python list would happily have held the float, which is why habits carried over from lists get people into trouble.

If you need decimals, declare that intent when the array is built—np.array([10, 20, 30], dtype=float)—or write 10.0 somewhere in the literal and let inference do it for you.