Numpy¶

$\newcommand{\re}{\mathbb{R}}$ Numpy provides many useful types and methods to perform numerical computations including vectors, matrices and linear algebra.

First import the numpy module.

You only need to do this once, usually at the beginning of your code.

To access methods, etc. inside this module, use dot operator.

It is common practice to import numpy like this.

Note that np acts as an alias or short-hand for numpy.

1-d arrays¶

Create an array of zeros

These one dimensional arrays are of type ndarray

Create an array of ones

Create and add two arrays

Get the size of array

Get the shape of an array

We see that these are arrays of reals by default. We can specify the type

linspace¶

This behaves same way as Matlab's linspace function.

Generate 10 uniformly spaced numbers in [1,10]

Note that this includes the end points 1 and 10. The output of linspace is an ndarray of floats.

x = linspace(a,b,n) is such that x is an array of n elements

x[i] = a + i*h,   i=0,1,2,...,n-1,    h = (b-a)/(n-1)

so that

x[0] = a,   x[n-1] = x[-1] = b

Note that we get last element of array using x[-1]; similarly, the last but one is given by x[-2], etc.

arange¶

For $m < n$, arange(m,n) returns

$$ m, m+1, \ldots, n-1 $$

Note the last value $n$ is not included.

We can specify a step size

If the arguments are float, the returned value is also float.

Beware of pitfalls - 1¶

Create an array of ones.

Maybe we want set all elements to zero, so we might try this

x has changed from an array to a scalar !!! In Python, the type of a variable is determined by what we assign to it. The correct way to set zeroes is this.

Beware of pitfalls - 2¶

Why did y change ? This happened because when we do

y = x

then y is just a pointer to x, so that changing x changes y also. If we want y to be an independent copy of x then do this

2-d arrays¶

2-d arrays can be considered as matrices, though Numpy has a separate matrix class.

Create an array of zeros

The size must be specified as a tuple or list.

We can access the elements by two indices A[i,j]: here i is the row index and j is the column index.

 -----------j-------->
 |  [[0. 0. 0. 0. 0.]
 |   [0. 0. 0. 0. 0.]
 i   [0. 0. 0. 0. 0.]
 |   [0. 0. 0. 0. 0.]
 v   [0. 0. 0. 0. 0.]]

The top, left element is A[0,0].

We can modify the elements of an array.

Create an array of ones

Create identity matrix

Create an array by specifying its elements

Create a random array and inspect its shape

To print the elements of an array, we need two for loops, one over rows and another over the columns.

To transpose a 2-d array

Diagonal matrix¶

Create

$$ A = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 2 & 0 & 0 \\ 0 & 0 & 3 & 0 \\ 0 & 0 & 0 & 4 \end{bmatrix} $$

Create $n \times n$ tri-diagonal matrix

empty_like, etc.¶

Sometimes we have an existing array and we want to create a new array of the same shape and type, but without initializing its elements.

1-D array is neither row nor column vector¶

Create a row vector

Create a column vector

Functions like zeros and ones can return row/column vector rather than array.

Row/column vectors are like row/column matrix. We have to use two indices to access the elements of such row/column vectors.

Here we access elements of a column vector.

x[0] gives an array for the first row.

Accessing portions of arrays¶

Array of 10 elements

x[0]	x[1]	x[2]	x[3]	x[4]	x[5]	x[6]	x[7]	x[8]	x[9]
0	1	2	3	4	5	6	7	8	9
x[-10]	x[-9]	x[-8]	x[-7]	x[-6]	x[-5]	x[-4]	x[-3]	x[-2]	x[-1]

Get elements x[2],...,x[5]

Hence x[m:n] gives the elements x[m],x[m+1],...,x[n-1].

Get elements from x[5] upto the last element

Get the last element

Get element x[5] upto last but one element

Access every alternate element of array starting from beginning x[0]

and starting from x[1]

These operations work on multi dimensional arrays also.

Arithmetic operations on arrays¶

Some arithmetic operations act element-wise

If A and x are arrays, then A*x does not give matrix-vector product. For that use dot

or equivalently

In Python3 versions, we can use @ to achieve matrix operations

We can of course do matrix-matrix products using dot or @

On vectors, @ performs dot product, i.e., x@y = dot(x,y)

Some matrix/vector functions¶

Note that np.sum(x*x) gives

$$ \sum_{i=0}^{n-1} x_i^2 $$

We can compute vector norms using numpy.linalg.norm (also see scipy.linalg.norm)

or import the linalg module and use methods inside it.

These methods work on 2-d arrays treated as matrices.

Example: Matrix-vector product¶

Let

$$ x \in \re^n, \qquad A \in \re^{m \times n}, \qquad y = A x \in \re^m $$

We can compute this element-wise

$$ y_i = \sum_{j=0}^{n-1} A_{ij} x_j, \qquad 0 \le i \le m-1 $$

We can verify that our result is correct by computing $\| y - A x\|_\infty$

We can also compute the product column-wise. Let

$$ A_{:,j} = \textrm{j'th column of A} \in \re^m $$

Then the matrix-vector product can also be written as

$$ y = \sum_{j=0}^{n-1} A_{:,j} x_j $$

Warning: This may have inefficient memory access since by default, numpy arrays have column-major ordering, see below.

Example: Matrix-Matrix product¶

If $A \in \re^{m\times n}$ and $B \in \re^{n \times p}$ then $C = AB \in \re^{m \times p}$ is given by

$$ C_{ij} = \sum_{k=0}^{n-1} A_{ik} B_{kj}, \qquad 0 \le i \le m-1, \quad 0 \le j \le p-1 $$

Let us verify the result is correct by computing

$$ |C-AB|_\infty = \max_{i,j} |C_{ij} - (AB)_{ij}| $$

using @ operator

Another view-point is the following

$$ C_{ij} = (\textrm{i'th row of A}) \cdot (\textrm{j'th column of B}) $$

Outer product of two vectors¶

Given $a \in \re^m$ and $b \in \re^n$, both considered as column vectors, their outer product

$$ A = a b^\top = \begin{bmatrix} | \\ a \\ | \end{bmatrix} \begin{bmatrix} - b^\top - \end{bmatrix} = \begin{bmatrix} | & | & & | \\ a b_0 & a b_1 & \ldots & a b_{n-1} \\ | & | & & | \end{bmatrix} \in \re^{m \times n} $$

is an $m \times n$ matrix with elements

$$ A_{ij} = a_i b_j, \qquad 0 \le i \le m-1, \quad 0 \le j \le n-1 $$

Math functions¶

Numpy provides standard functions like sin, cos, log, etc. which can act on arrays in an element-by-element manner. This is not the case for functions in math module, which can only take scalar arguments.

Memory ordering in arrays*¶

By default, the ordering is same as in C/C++, the inner-most index is the fastest running one. For example, if we have an array of size (2,3), they are stored in memory in this order

a[0,0], a[0,1], a[0,2], a[1,0], a[1,1], a[1,2]

i.e.,

a[0,0] --> a[0,1] --> a[0,2]
                         |
   |----------<----------|
   |
a[1,0] --> a[1,1] --> a[1,2]

To get fortran style ordering, where the outer-most index is the fastest running one, which corresponds to the following layout in memory

a[0,0], a[1,0], a[0,1], a[1,1], a[0,2], a[1,2]

i.e.,

a[0,0]    -- a[0,1]    -- a[0,2]
   |     |     |      |      |  
   v     ^     v      ^      v
   |     |     |      |      |
a[1,0] --    a[1,1] --    a[1,2]

create like this

Tensor product array: meshgrid¶

If $x \in \re^m$ and $y \in \re^n$, then the output of meshgrid(x,y) is arranged like this

$$ X[i,j] = x[j], \qquad Y[i,j] = y[i], \qquad 0 \le i \le n-1, \quad 0 \le j \le m-1 $$

If we want the following arrangement

$$ X[i,j] = x[i], \qquad Y[i,j] = y[j], \qquad 0 \le i \le m-1, \quad 0 \le j \le n-1 $$

we have to do the following

The second form is useful when working with finite difference schemes on Cartesian grids, where we want to use i index running along x-axis and j index running along y-axis.

Reshaping arrays¶

Writing and reading files¶

Write an array to file

Check the contents of the file in your terminal

cat A.txt

We can also print it from within the notebook

Write two 1-D arrays as columns into file

We can control the number of decimals saved, and use scientific notation

We can read an existing file like this

Measuring time¶

Loops are slow in Python and it is good to avoid them. The following example shows how to time code execution.

Create some random matrices

First we add two matrices with an explicit for loop.

Now we use the + operator.

The loop version is significantly slower.