Numpy arrays are more efficient than lists because all elements of numpy arrays are of the same pre-determined type. Numpy also provides efficient ufuncs (universal functions) which are vectorized functions that loop over array elements with loops pre-compiled in C. Numpy also exhibits some syntactic features that a mathematician may consider a nuisance, but knowing how to work with numpy is key for almost all scientific computation in python. It takes some practice to be comfortable with numpy and practice is what we are aiming for in this activity. This activity is structured as a list of questions. You will be in a position to appreciate these questions (and their answers) after going through this lecture's prerequistes on numpy yourself.

In [1]:

```
import numpy as np
import math
```

In [2]:

```
A = [0.1, 1.3, 0.4, 0.5] # list
a = np.array(A) # numpy array
type(a), type(A)
```

Out[2]:

Here is how you find out the common data type of elements of a numpy array (and there is no such analogue for list, since list elements can be of different types).

In [3]:

```
a.dtype # a's common element type (A.dtype is undefined!)
```

Out[3]:

`2*a`

and `2*A`

?¶In [4]:

```
2*a
```

Out[4]:

In [5]:

```
2*A
```

Out[5]:

Here is one option:

```
[sin(x[i]) * exp(-x[i]) for i in range(n)]
```

And here is another:

```
np.sin(x) * np.exp(-x)
```

Which is better?

In [6]:

```
n = 100000
x = np.linspace(0, 2*np.pi, n)
```

In [7]:

```
# list comprehension
%timeit y = [math.sin(x[i]) * math.exp(-x[i]) for i in range(n)]
# use numpy ufuncs
%timeit y = np.sin(x) * np.exp(-x)
```

The functions `np.sin`

and `np.exp`

are examples of numpy's universal functions (ufuncs) that act directly on arrays. While `np.sin`

are unary ufunc, there are many binary ufuncs like `np.add`

: when you write `x+y`

, you are actually calling the binary ufunc `np.add(x, y)`

. Ufuncs are *vectorized* functions.

Vectorization refers to one or both of the following, depending on context:

A convenience feature: Apply an operation to all elements of a collection at once.

A performance feature: Use hardware instruction sets to execute single instruction on multiple data (SIMD).

A numpy ufunc like `np.sin`

is vectorized in both the above senses: (1) you can apply it directly to an array thus avoiding python loops (unlike `math.sin`

which can be applied to just a single value), and (2) numpy's C implementation of `np.sin`

uses some SIMD instruction sets, allowing you to get automatic speed up when running on hardware that supports the instructions. If the latter still sounds mysterious, here is more explanation than you probably need: most chips today come with a SIMD instruction to process 4 numbers (float64) at once (and fancier chips can do more), so a loop over an array of N floats can finish in N/4 iterations if you utilize that instruction.

To examine the difference between the convenience feature and the performance feature, consider the following function, which is written to apply to just one number:

In [8]:

```
def f(v): # apply f to one scalar value v
return math.sin(v) * math.exp(-v)
```

To gain (1), the convenience feature, there are at least two options other than using ufuncs:

*a) Use map*

A function `f`

acting on a scalar value can be made into a function that acts on a vector of values using the functional programming tool `map(f, x)`

, which
returns an iterator that applies f to every element of x.

In [9]:

```
vectorizedf = map(f, x) # apply same f to a vector of values x
```

*b) Use numpy's vectorize*

In [10]:

```
F = np.vectorize(f) # F can be applied to a array x
```

Both options (a) and (b) provide the convenience feature (1), letting you avoid python loops, and allows you to write expressive short codes.

However, neither option (a) nor (b) gives you the full performance of numpy's ufuncs.

In [11]:

```
# use map
%timeit y = list(map(f, x))
# use numpy's vectorize
%timeit y = F(x)
# use numpy's ufunc
%timeit y = np.sin(x) * np.exp(-x)
```

`range`

as efficient as `np.arange`

?¶In [12]:

```
%timeit for x in range(1000000): x**3
%timeit for x in np.arange(1000000): x**3
```

There was a time (in Python 2) when `range`

was not as efficient, but those times have passed.

A slice

`a[b:e:s]`

of`a`

refers to the array of elements from the**b**eginning index`b`

(included) till**e**nding index`e`

(excluded),**s**tepping`s`

elements.The defaults

`b=0`

,`e=len(a)`

, and`s=1`

may be omitted in a slice specification.Negative indices count from the end of the array:

`a[-1]`

is the last element of`a`

and`a[-k] = a[len(a)-k]`

.Positive indices count from the begining as usual.

In [13]:

```
a = np.random.randint(0,9,5)
a
```

Out[13]:

If you have understood these, then you should be able to say what the expected results are from each of the following statements.

```
a[::]
a[-1]
a[len(a)-1]
a[-3:]
a[-4:-1:2]
slice = range(-4,-1,2)
a[-4:-1:2], a[slice]
```

Verify your answers:

In [14]:

```
a[::]
```

Out[14]:

In [15]:

```
a[-3:]
```

Out[15]:

In [16]:

```
a[-1], a[len(a)-1]
```

Out[16]:

In [17]:

```
a[-4:-1:2]
```

Out[17]:

In [18]:

```
slice = range(-4,-1,2) # Think of b:e:s specification as a range.
a[-4:-1:2], a[slice] # In older versions, a[slice] may not work
# but will work with slice=arange(-4,-1,2).
```

Out[18]:

`=`

does?¶In [19]:

```
a = np.array([1,2,3])
b = np.array([3,4,5,6])
```

After assigning `a`

to `b`

by `=`

, what happens when you change an element of `a`

?

In [20]:

```
a = b
a[0] = 1
a
```

Out[20]:

We certainly expected the 3 to become 1 in `a`

. Did you also expect the following?

In [21]:

```
b
```

Out[21]:

If this surprises you, pay close attention to the next question.

In most languages, each variable has its own memory address. For example consider this simple C++ code (ignore it if you don't know C++).

```
#include <vector>
std::vector<int> a{1,2,3}, b{3,4,5,6};
// Objects a and b each have their own memory addresses.
// Assignment a=b copies contents of b's memory into a.
a = b;
// a's memory address has not changed, but its contents have.
```

If you have programmed in C or C++, you might have gotten used to variables being permanently linked to their memory locations.

In contrast, python variables are just names. In python, variables like `a`

and `b`

are *names* which are not associated to fixed memory addresses. Names can be *bound* to one object in memory, and later to another. Multiple names can be bound to the same object (sometimes known as *aliasing* in other languages). The upshot of this is that in python, the assignment "`=`

" changes names, but need not copy memory contents.

In [22]:

```
a = np.array([1,2,3]) # This is Object1 and "a" is a name for it.
b = np.array([3,4,5,6]) # This is Object2 and "b" is a name for it.
```

We can double check that these objects occupy different memory locations using python's `id`

function.

In [23]:

```
id(a), id(b)
```

Out[23]:

Consider what happens when you say `a=b`

:

In [24]:

```
a = b # a is no longer a name for Object1, it is now a name for Object2.
```

In [25]:

```
id(a), id(b)
```

Out[25]:

Names `a`

and `b`

are now both bound to the same "Object2". (And we no longer have a way to access "Object1"!) Now, when we change `a`

, or when we change `b`

, we are really changing the same object.

In [26]:

```
a = np.array([1,2,3]) # Object1
b = np.array([3,4,5,6]) # Object2
a = b.copy() # Copies Object2, and binds a to the copy
a[0] = 2 # Only the copied (new) object is changed
```

In [27]:

```
a, b
```

Out[27]:

Of course, numpy is all about vectors and matrices (and even higher-order tensors). Two-dimensional data, or tabular data, or matrix data of the form $$ \left\lbrack\begin{array}{ccc} A_{0,0} & \cdots & A_{0,n-1}\\ \vdots & \ddots & \vdots\\ A_{m-1,0} & \cdots & A_{m-1,n-1} \end{array}\right\rbrack $$ can be represented in python

- either as list of lists
- or as a numpy array.

The numpy `array`

is more efficient than list of lists and has constructs for some matrix operations. (Note that you might find a `numpy.matrix`

class, distinct from the `array`

class, in some older codes, but be warned that it is deprecated. Due to problems arising from mixing `matrix`

and `array`

objects in python codes, we will *not* use the deprecated `matrix`

class in this course. You should not use it in work you turn in.)

In [28]:

```
Amat = [[1,2],
[3,4]]
Amat
```

Out[28]:

In [29]:

```
amat = np.array(Amat)
amat
```

Out[29]:

In [30]:

```
type(A), type(a)
```

Out[30]:

Note that 2D and 1D numpy arrays are of the same type called `numpy.ndarray`

.

What is the difference between `2*Amat`

and `2*amat`

, for the objects `Amat`

(list of lists) and `amat`

(numpy array) just made above?

In [31]:

```
2*Amat
```

Out[31]:

In [32]:

```
2*amat
```

Out[32]:

In [33]:

```
amat
```

Out[33]:

In [34]:

```
amat * amat
```

Out[34]:

Look at the output: *is this really matrix multiplication?!* This is one thing that drives mathematicians crazy when they look at `numpy`

for the first time. Mathematicians want the multiplication operator `*`

to mean matrix multiplication when it is applied to numpy arrays. Unfortunately python's default `*`

does **element-by-element multiplication**, not matrix multiplication. Since the first proposal for numpy, decades ago, many battles have been waged to resolve this embarrassment.

Finally, a few years ago, there came some good news. Since Python 3.5, the `@`

symbol was dedicated to mean the matrix multiplication operator. You can read more about it at PEP 465.

In [35]:

```
import sys
print(sys.version) # check if you have version >= 3.5 before trying @
```

In [36]:

```
amat @ amat
```

Out[36]:

Naturally, many of us needed matrix multiplication before the `@`

came along, so as you must have guessed, there is another way to do matrix multiplication:

In [37]:

```
np.dot(amat, amat) # dot(A,B) = matrix A multiplied by matrix B
```

Out[37]:

In [38]:

```
amat.dot(amat)
```

Out[38]:

I think you will agree with me that this is not as neat as `@`

.

You should know that the embarrassment continues in matrix powers. If you thought `amat ** 2`

should give you a matrix power equaling the product of `amat`

with itself, think again.

In [39]:

```
amat**2 # not equal to matrix power !!
```

Out[39]:

Numpy provides a `matrix_power`

routine to compute $M^n$ for matrices $M$ (and integers $n$).

In [40]:

```
np.linalg.matrix_power(amat, 2)
```

Out[40]:

It does the job, but leaves elegance by the wayside.

Slicing in two-dimensional arrays is similar to slicing one-dimensional arrays.
If `rslice`

and `cslice`

are 1D slices (or ranges) like the ones we used for one-dimensional arrays, then when applied to a 2D array `A`

,

```
A[rslice, cslice]
```

the result is a submatrix of `A`

using row indices in `rslice`

and column indices in `cslice`

.

In [41]:

```
A = np.array([[7, 8, 5, 1], [2, 5, 5, 2], [9, 6, 8, 9]])
A
```

Out[41]:

In [42]:

```
A[1, :], A[:, 2]
```

Out[42]:

In [43]:

```
A[:3:2, :3]
```

Out[43]:

Like other programming facilities, numpy stores 2D array data internally as a 1D array, in order to get a contiguous memory block for convenient storage. For 2D arrays, Fortran and Matlab uses column-major ordering, while C uses row-major ordering. For example, the 2D array

```
[[7, 8, 5, 1],
[2, 5, 5, 2],
[9, 6, 8, 9]]
```

in *row-major ordering* looks like

`7, 8, 5, 1, 2, 5, 5, 2, 9, 6, 8, 9`

while in *column-major ordering*, it looks as follows.

`7, 2, 9, 8, 5, 6, 5, 5, 8, 1, 2, 9`

Numpy, by default, stores arrays in row-major ordering (like C). This thinking is reflected in some numpy's methods: e.g., when you ask numpy to reshape or flatten a array, the result is what you expect as if it were stored in row-major ordering.

In [44]:

```
M = np.array([[7, 8, 5, 1], [2, 5, 5, 2], [9, 6, 8, 9]])
M
```

Out[44]:

In [45]:

```
M.reshape(2, 6) # Just a different view of the same data
```

Out[45]:

In [46]:

```
M.ravel() # The 1D data of M in row-major ordering
```

Out[46]:

But the actual situation is more complicated since numpy allows users to *override the default storage* ordering. You can decide to store an array like in C or like in Fortran. Here is how to store the same array in Fortran's column-major ordering.

In [47]:

```
A = np.array(M, order='F')
A
```

Out[47]:

Obviously, it is the same matrix. How the data is stored internally is mostly immaterial (except for some performance optimizations). The behavior (of most) of numpy methods does not change even if the user has opted to store the array in a different ordering. If you really need to see a matrix's internal ordering, you can do so by calling the `ravel`

method with keyword argument `order='A'`

.

In [48]:

```
A.ravel(order='A') # A's internal ordering is Fortran style
```

Out[48]:

In [49]:

```
M.ravel(order='A') # M's internal ordering is default C-style
```

Out[49]:

If a numpy array is given indices that are boolean (instead of integers), then rows or columns are selected based on `True`

indices. This is called **masking**. It is very useful together with vectorized conditionals.

In [50]:

```
N = np.arange(25).reshape(5,5)
N
```

Out[50]:

How will you isolate elements in `N`

whose value is between 7 and 18?

In [51]:

```
mask = (N>7) & (N<18)
mask
```

Out[51]:

These are the elements we needed:

In [52]:

```
N[mask]
```

Out[52]:

And these are their locations:

In [53]:

```
i, j = np.where(mask) # Returns i and j indices where mask[i,j] is True.
i, j # 1st True value of mask is at i[0],j[0],
# 2nd True value of mask is at i[1],j[1], etc.
```

Out[53]:

Numpy can work with general-dimensional arrays, not just 1D or 2D arrays. For an `n`

-dimensional array, the `shape`

of a numpy array is a tuple of `n`

integers giving the sizes in each dimension.

In [54]:

```
data = np.random.randint(low=0, high=10, size=30) # 1D array
```

In [55]:

```
T2 = np.reshape(data, (6, 5)) # 2D array
T2
```

Out[55]:

In [56]:

```
T3 = np.reshape(data, (2, 3, 5)) # 3D array
T3
```

Out[56]:

In [57]:

```
print('T3 is a ', T3.ndim, 'dimensional array of shape ', T3.shape)
print('T2 is a ', T2.ndim, 'dimensional array of shape ', T2.shape)
print('data is a ', data.ndim, 'dimensional array of shape ', data.shape)
```

Here are a few other features of numpy arrays to note:

- Every numpy array has attributes
`ndim`

and`shape`

. - A scalar
`c`

, or`np.array(c)`

is considered to have`ndim=0`

and`shape=()`

. - A vector of length
`n`

, when viewed as a row vector has`ndim=1`

and`shape=(n,)`

. - A vector of length
`n`

, when viewed as a column vector has`ndim=2`

and`shape=(n, 1)`

. - You can convert a row vector
`a`

to a column vector by`a[:, np.newaxis]`

. - Use
`newaxis`

to add a new dimension, e.g.,`T3[:, :, np.newaxis, :]`

has`shape=(2, 3, 1, 5)`

.

In mathematics, it would be an illegal operation to add matrices of different shapes. But it is not surprising that we would want to: e.g., viewing the number 10 as a 1x1 matrix and considering a matrix `A`

of any other size, wouldn't it be nice to say `10 + A`

to add 10 to all elements of A? Wouldn't it also be nice to be able to use `+`

to add a vector to all columns of a matrix with more than one columns? All this and more is made possible in numpy by **broadcasting** rules, which extend the possibilities of vectorized operations. A very clear explanation of broadcasting is in [JV-H].

To see if you can add up (or apply another binary ufunc) differently shaped arrays, follow this algorithm, which uses the `ndim`

and `shape`

attributes we just saw.

Step 1. If two arrays differ in their `ndim`

, then revise the `shape`

of the one with lower `ndim`

by prepending 1 on the left until the `ndim`

s are equal.

Step 2. If `shape[i]`

(after its possible revision from Step 1) of the two arrays are unequal and one of them equals 1, then increase the latter to match the other `shape[i]`

.

If the resulting revised `shape`

s of the arrays are still unequal, then broadcasting fails and you can't add them. In Step 1, when we increase `ndim`

by prepending a 1 to `shape`

, we are not really changing the array: we are just imagining it with one extra dimension. In Step 2, when we increase `shape[i]`

from 1 to something larger, we are imagining the array elements repeated along the `i`

-th dimension (without actually performing an operation copying the elements). Here are a few examples to illustrate these rules.

*Example 1:*

```
a + b = [1, 8, 3] + [1]
[0, 6, 5] [8]
a.ndim=2 b.ndim=2 <= No need for Step 1 modification
a.shape=(2, 3) b.shape=(2, 1) <= Apply Step 2 to equalize shape
Step 2: a.shape=(2, 3) b.shape=(2, 3)
a + b = [1, 8, 3] + [1, 1, 1] = [ 2, 9, 4]
[0, 6, 5] [8, 8, 8] [ 8, 14, 13]
```

*Example 2:*

```
a + b = [1, 8, 3] + [1]
[0, 6, 5]
a.ndim=2 b.ndim=0 <= Apply Step 1 to equalize ndim
a.shape=(2, 3) b.shape=() (prepend 1 until ndim equalizes)
Step 1: a.shape=(2, 3) b.shape=(1, 1) <= Next apply Step 2 to equalize shape
Step 2: a.shape=(2, 3) b.shape=(2, 3)
a + b = [1, 8, 3] + [1, 1, 1] = [ 2, 9, 4]
[0, 6, 5] [1, 1, 1] [ 1, 7, 5]
```

*Example 3:*

```
a + b = [1, 8, 3] + [1, 3]
[0, 6, 5]
a.ndim=2 b.ndim=1 <= Apply Step 1 to equalize ndim
a.shape=(2, 3) b.shape=(2, )
Step 1: a.shape=(2, 3) b.shape=(1, 2) <= Next apply Step 2 to equalize shape
Step 2: a.shape=(2, 3) b.shape=(2, 2) <= Still unequal: broadcasting fails
```

As a simple exercise to further fix ideas, follow the above procedure and see if you can explain whether broadcasting rules apply, or not, to the following (with `T2`

and `T3`

as set previously).

`T2 + T3`

`T3 + 1`

`T2[:3,:] + T3`

Author: Jay Gopalakrishnan

License: ©2020. CC-BY-SA

$\ll$Table of Contents