Lesson A11 – NumPy

This tutorial will focus on the NumPy module. NumPy is fundamental for scientific computing and will provide the basis for our everyday work, as well as offer advanced functionalities for numerical operations.

Refer also to the NumPy Quickstart tutorial online or to

“Travis E. Oliphant. 2015. Guide to NumPy (2nd. ed.). CreateSpace Independent Publishing Platform, North Charleston, SC, USA.”

if you prefer books.

To import the module run the cell below. NumPy is a third-party Python module that you might need to install first (e.g. using the conda package manager with conda install numpy).

import numpy as np

We can now call any function x implemented into NumPy with np.x(). Lets generate some lists to get started. You should aready know about lists. Recall that a list in Python is a mutable, ordered collection of objects with no restriction on the type.

list1 = [1, 3, 5]
list2 = [0.0, 2.0, 4.0]

For numerical operations, this data type is of limited use.

Arrays basics

We want to turn these lists into NumPy arrays. Arrays are a type of collection better suited for certain tasks and are essential for our day-to-day-work. They offer a lot of functionality. To convert these lists to arrays we run:

A = np.asarray(list1)
B = np.asarray(list2)
print('A = ', A, type(A))
print('B = ', B, type(B))

A =  [1 3 5] <class 'numpy.ndarray'>
B =  [0. 2. 4.] <class 'numpy.ndarray'>

We could have also created a NumPy array directly.

A = np.array([1, 3, 5])
print('A = ', A, type(A))

A =  [1 3 5] <class 'numpy.ndarray'>

On a first glance, the arrays we created look like lists when printed, but they are in fact very different. The type of objects collected in an array is constrained to one specific type. We can discover the data type associated with an array from the dtype array-attribute.

print('Type stored in A: ', A, A.dtype)
print('Type stored in B: ', B, B.dtype)

Type stored in A:  [1 3 5] int64
Type stored in B:  [0. 2. 4.] float64

Note: Arrays are mutable, ordered sequences in which all stored objects need to have the same type.

Advanced: The Python standard library also offers a type of collection called array. These arrays are similar to the arrays from the NumPy module, with respect to the data type restriction, but they are handled differently. You can use the array module with import array. Do not confuse arrays of the type numpy.ndarray with those of the type array.array.

Because of the restriction to only one type, lists and other sequences can only be converted to NumPy arrays straight away if they fulfill this criterion. Otherwise the objects they contain are transformed into a common basic type if possible.

flexible_collection = [2, True, "x"]
np.asarray(flexible_collection)  # Converted to string object

array(['2', 'True', 'x'], dtype='<U21')

flexible_collection = [2, True, "x", None]
np.asarray(flexible_collection)  # Converted to `object`

array([2, True, 'x', None], dtype=object)

To quite a large extend arrays can be treated just as other sequences. We can for example iterate over them.

A = np.array([1, 3, 5])
for i in A:
    print(i**2)

1
9
25

Or we can access their elements by indexing.

print(A)        # print the whole string
print(A[0])     # print the first character of the string
print(A[-1])    # print the last character of the string
print(A[:2])    # print the first two characters of the string
print(A[1:-1])  # print the second to the last character of the string
print(A[::2])   # print every second character
print(A[::-1])  # print every character in reversed order
print(A[10])    # raises IndexError

[1 3 5]
1
5
[1 3]
[3]
[1 5]
[5 3 1]

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-9-48fc44548521> in <module>
      6 print(A[::2])   # print every second character
      7 print(A[::-1])  # print every character in reversed order
----> 8 print(A[10])    # raises IndexError

IndexError: index 10 is out of bounds for axis 0 with size 3

The indexing capabilities of arrays go, however, even further. We can for example use another sequence (e.g. another array or list) to specifiy the indices of elements we want to get.

# Using a sequence of indices to index an array
A = np.array([1, 3, 5, 7])  # Array of elements
B = [1, 3]                  # List of indices
print(A[B]) # Get all elements in A with indices given in B

[3 7]

This is not working for lists for example.

# Trying to use a sequence of indices to index a list
A = [1, 3, 5, 7]  # List of elements
B = [1, 3]        # List of indices
print(A[B])  # Can't get elements in A with indices given in B

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-11-935edb6c284e> in <module>
      2 A = [1, 3, 5, 7]  # List of elements
      3 B = [1, 3]        # List of indices
----> 4 print(A[B])  # Can't get elements in A with indices given in B

TypeError: list indices must be integers or slices, not list

Also very useful, we can use a sequence of objects of the boolean type (True and False) with the same length as the array we want to index to get only specific values.

# Using a boolean sequence of indices to index an array
A = np.asarray([1, 3, 5, 7])
B = [True, True, False, True]
# A and B must be of same length
print(A[B])

[1 3 7]

# Trying to use a boolean sequence of indices to index a list
A = [1, 3, 5, 7]
B = [True, True, False, True]
print(A[B])

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-13-5b4bb0a00302> in <module>
      2 A = [1, 3, 5, 7]
      3 B = [True, True, False, True]
----> 4 print(A[B])

TypeError: list indices must be integers or slices, not list

This extends to testing an array against a condition to get a boolean sequence for indexing (see also Boolean expressions and np.where() to find them).

A = np.asarray([1, 3, 5, 7])
print(A[A > 1])

[3 5 7]

One of the most striking advantages of arrays over other sequences, is that we can perform certain operations on all of their elements at the same time. We can add, subtract, multiply or divide scalars and the arrays. The operation will be performed element-wise.

A = np.asarray([1, 3, 5, 7])
print('A + 5.0 = ', A + 5.0)
print('A - 5,0 = ', A - 5.0)
print('A * 5.0 = ', A * 5.0)
print('A / 5.0 = ', A / 5.0)

A + 5.0 =  [ 6.  8. 10. 12.]
A - 5,0 =  [-4. -2.  0.  2.]
A * 5.0 =  [ 5. 15. 25. 35.]
A / 5.0 =  [0.2 0.6 1.  1.4]

We can also add, subtract, multiply or divide arrays (of same length/size). The operation will be performed again element-wise.

A = np.asarray([1, 3, 5, 7])
B = np.asarray([0, 2.0, 4.0, 8.5])
print('A + B = ', A + B)
print('A - B = ', A - B)
print('A * B = ', A * B)
print('A / B = ', A / B)

A + B =  [ 1.   5.   9.  15.5]
A - B =  [ 1.   1.   1.  -1.5]
A * B =  [ 0.   6.  20.  59.5]
A / B =  [       inf 1.5        1.25       0.82352941]

<ipython-input-16-42687086e92f>:6: RuntimeWarning: divide by zero encountered in true_divide
  print('A / B = ', A / B)

The division by zero in the last line is the reason behind the warning message (RuntimeWarning). Note that a division by zero outside NumPy would result in a ZeroDivisionError. NumPy can handle this and has a type to represent infinity (np.inf).

1 / 0

---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)
<ipython-input-17-bc757c3fda29> in <module>
----> 1 1 / 0

ZeroDivisionError: division by zero

The operations above can be translated to vector operations in linear algebra. Let’s define \(\mathbf{a}\) and \(\mathbf{b}\) as the vectors corresponding to the NumPy arrays A and B.

\begin{equation*} \mathbf{a} = \begin{pmatrix} 1\\ 3\\ 5\\ 7 \end{pmatrix},~~~~~~~ \mathbf{b} = \begin{pmatrix} 0\\ 2.0\\ 4.0\\ 8.57 \end{pmatrix} \end{equation*}

Let \(\mathbf{e} = \mathbf{1}\) be a vector of all ones and \(\mathbf{I}\) the identity matrix.

\begin{equation*} \mathbf{e} = \begin{pmatrix} 1\\ 1\\ 1\\ 1 \end{pmatrix},~~~~~~~ \mathbf{I} = \begin{pmatrix} 1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1 \end{pmatrix} \end{equation*}

Then, the NumPy operations have the following algebraic expressions:

A + B \(\rightarrow\) \(\mathbf{a} + \mathbf{b}\)

A - B \(\rightarrow\) \(\mathbf{a} - \mathbf{b}\)

A * 5 \(\rightarrow\) \(5 \mathbf{a}\)

A / 5 \(\rightarrow\) \(\frac{1}{5} \mathbf{a}\)

A + 5 \(\rightarrow\) \(\mathbf{a} + 5\cdot\mathbf{e}\)

A - 5 \(\rightarrow\) \(\mathbf{a} - 5\cdot\mathbf{e}\)

A * B \(\rightarrow\) \(\mathbf{a} \cdot \mathbf{I} \cdot \mathbf{b}\)

A / B \(\rightarrow\) \(\mathbf{a} \cdot (\mathbf{I} \cdot \mathbf{b})^{-1}\)

With \((\mathbf{I} \cdot \mathbf{b})^{-1}\) being the inverse of \(\mathbf{I} \cdot \mathbf{b}\).

Note: Be careful, A * B is not the scalar product \(\mathbf{a} \cdot \mathbf{b}\). This you get with np.dot(A, B) or A @ B.

Arrays come with many handy methods like sum(), mean(), std() etc. Here are some examples.

Note: Remember that you can get an overview over all these methods by using dir(object). You can also explore them by typing "A." and then pressing tab to use the interactive help in a Jupyter notebook.

print("A.sum() = ", A.sum())
# The largest Entry of the Array
print("A.max() = ", A.max())
# The index of the largest entry
print("A.argmax() = ", A.argmax())

A.sum() =  16
A.max() =  7
A.argmax() =  3

Arrays can be multidimensional. This is how you represent matrices in NumPy. The dimensionality of an array can be viewed with the shape property. We are going to make a 2D-array C, which will be the array made from a list of lists.

list1 = [1, 3, 5]
list2 = [0.0, 2.0, 4.0]

C = np.array([
    list1,
    list2,
    ])

print("C = ")
print(C)
print("\n")
print('C.shape (rows, columns) = ', C.shape)
print("\n")
print('C[0] = ', C[0])
print('C[1] = ', C[1])
print('C[:, 0] = ', C[:, 0])  # [rowindex, columnindex]

C =
[[1. 3. 5.]
 [0. 2. 4.]]


C.shape (rows, columns) =  (2, 3)


C[0] =  [1. 3. 5.]
C[1] =  [0. 2. 4.]
C[:, 0] =  [1. 0.]

Note that this is fundamentally different to having a nested list of lists.

nested_list = [
    list1,
    list2,
    ]

print("nested_list[0][0] (nested index) =>", nested_list[0][0])
nested_list[:, 0]  # raises TypeError

nested_list[0][0] (nested index) => 1

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-20-e19173c88ae2> in <module>
      5
      6 print("nested_list[0][0] (nested index) =>", nested_list[0][0])
----> 7 nested_list[:, 0]  # raises TypeError

TypeError: list indices must be integers or slices, not tuple

If you depend on fast row- or column-wise indexing, e.g. when working with matrices, NumPy arrays can make your life much easier.

Boolean expressions and np.where() to find them

Another essential function that deserves extra mentioning is np.where(). In the most simple case, it takes a boolean sequence, i.e. a sequence of True and False, as an input and returns the indices where the sequence is true.

help(np.where)

Help on function where in module numpy:

where(...)
    where(condition, [x, y])

    Return elements chosen from `x` or `y` depending on `condition`.

    .. note::
        When only `condition` is provided, this function is a shorthand for
        ``np.asarray(condition).nonzero()``. Using `nonzero` directly should be
        preferred, as it behaves correctly for subclasses. The rest of this
        documentation covers only the case where all three arguments are
        provided.

    Parameters
    ----------
    condition : array_like, bool
        Where True, yield `x`, otherwise yield `y`.
    x, y : array_like
        Values from which to choose. `x`, `y` and `condition` need to be
        broadcastable to some shape.

    Returns
    -------
    out : ndarray
        An array with elements from `x` where `condition` is True, and elements
        from `y` elsewhere.

    See Also
    --------
    choose
    nonzero : The function that is called when x and y are omitted

    Notes
    -----
    If all the arrays are 1-D, `where` is equivalent to::

        [xv if c else yv
         for c, xv, yv in zip(condition, x, y)]

    Examples
    --------
    >>> a = np.arange(10)
    >>> a
    array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
    >>> np.where(a < 5, a, 10*a)
    array([ 0,  1,  2,  3,  4, 50, 60, 70, 80, 90])

    This can be used on multidimensional arrays too:

    >>> np.where([[True, False], [True, True]],
    ...          [[1, 2], [3, 4]],
    ...          [[9, 8], [7, 6]])
    array([[1, 8],
           [3, 4]])

    The shapes of x, y, and the condition are broadcast together:

    >>> x, y = np.ogrid[:3, :4]
    >>> np.where(x < y, x, 10 + y)  # both x and 10+y are broadcast
    array([[10,  0,  0,  0],
           [10, 11,  1,  1],
           [10, 11, 12,  2]])

    >>> a = np.array([[0, 1, 2],
    ...               [0, 2, 4],
    ...               [0, 3, 6]])
    >>> np.where(a < 4, a, -1)  # -1 is broadcast
    array([[ 0,  1,  2],
           [ 0,  2, -1],
           [ 0,  3, -1]])

np.where([True, False, True])  # returns (x-indices, y-indices)

(array([0, 2]),)

It is most often used when the boolean sequence is derived from a boolean expression involving an array.

A = np.array([1, 3, 5, 7])
print(A == 3)
print(np.where(A == 3))  # returns (x-indices, y-indices)

[False  True False False]
(array([1]),)

It is of course also working in two dimensions. Lets see where C is greater than or equal to 3.

list1 = [1, 3, 5]
list2 = [0.0, 2.0, 4.0]

C = np.array([
    list1,
    list2,
    ])

print(C)

[[1. 3. 5.]
 [0. 2. 4.]]

This can be also expressed as the corresponding matrix:

\begin{equation*} \mathbf{C} = \begin{pmatrix} 1& 3& 5\\ 0& 2& 4 \end{pmatrix} \end{equation*}

print(np.where(C >= 3))  # returns (x-indices, y-indices)

(array([0, 0, 1]), array([1, 2, 2]))

The result are two arrays with three indices each. We get two arrays because C is two-dimensional. Each array has three entries because in total there are three hits. This means that the statement is True for the indices C[0, 1], C[0, 2] and C[1, 2]. Lets check:

# All greater than three?
print("C[0, 1] =", C[0, 1])
print("C[0, 2] =", C[0, 2])
print("C[1, 2] =", C[1, 2])

C[0, 1] = 3.0
C[0, 2] = 5.0
C[1, 2] = 4.0

Since indices in Python are zero-based but they are 1-based in terms of linear algebra, we need to shift the indices by one to get the corresponding matrix elements.

\(\mathbf{C}_{12} = 3\)

\(\mathbf{C}_{13} = 5\)

\(\mathbf{C}_{23} = 4\)

Advanced

Optionally, we can pass two additional sequences to np.where from which values are returned instead of just the true-indices.

print(np.where([True, False, True], ["yes", "yes", "yes"], ["no", "no", "no"]))

['yes' 'no' 'yes']

Or simply:

np.where([True, False, True], "yes", "no")

array(['yes', 'no', 'yes'], dtype='<U3')

Special array initialisation

We have seen how to create arrays from sequences and from scratch. NumPy provides a handful of convenience functions to construct arrays filled with all kinds of values. There is for example the function np.linspace, which takes in principle three arguments (start, stop, num) and returns a serial array.

print(np.linspace(0, 10, 11))

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

We have also the possibilty to create arrays full of zeros, ones, or another value.

print(np.zeros(10))
print(np.ones(10))
print(np.full(10, "x"))

[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
['x' 'x' 'x' 'x' 'x' 'x' 'x' 'x' 'x' 'x']

Meshgrids

Consider the following: imagine you want to compute some function \(f\) that depends on spherical coordinates \(r\), \(\theta\) and \(\phi\). If you want a resolution of 5 unit points in each dimension we might naively run the following.

res = 5  # set resolution

r = np.linspace(0, 1, res)
theta = np.linspace(0, np.pi, res)
phi = np.linspace(-np.pi, np.pi, res)

print('r = ', r)
print('theta = ', theta)
print('phi = ', phi)

r =  [0.   0.25 0.5  0.75 1.  ]
theta =  [0.         0.78539816 1.57079633 2.35619449 3.14159265]
phi =  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]

Note how linspace produces evenly spaced numbers in defined ranges.

However, when we want to calculate \(f\) with the arrays r, theta and phi, the problem is that we need all possible combinations of the three input coordinates to fully describe a sphere. The solution is called np.meshgrid and works with two or more input arrays.

R, T, P = np.meshgrid(r, theta, phi)
print("R = ")
print(R)

print("\n")
print("T = ")
print(T)

print("\n")
print("P = ")
print(P)

R =
[[[0.   0.   0.   0.   0.  ]
  [0.25 0.25 0.25 0.25 0.25]
  [0.5  0.5  0.5  0.5  0.5 ]
  [0.75 0.75 0.75 0.75 0.75]
  [1.   1.   1.   1.   1.  ]]

 [[0.   0.   0.   0.   0.  ]
  [0.25 0.25 0.25 0.25 0.25]
  [0.5  0.5  0.5  0.5  0.5 ]
  [0.75 0.75 0.75 0.75 0.75]
  [1.   1.   1.   1.   1.  ]]

 [[0.   0.   0.   0.   0.  ]
  [0.25 0.25 0.25 0.25 0.25]
  [0.5  0.5  0.5  0.5  0.5 ]
  [0.75 0.75 0.75 0.75 0.75]
  [1.   1.   1.   1.   1.  ]]

 [[0.   0.   0.   0.   0.  ]
  [0.25 0.25 0.25 0.25 0.25]
  [0.5  0.5  0.5  0.5  0.5 ]
  [0.75 0.75 0.75 0.75 0.75]
  [1.   1.   1.   1.   1.  ]]

 [[0.   0.   0.   0.   0.  ]
  [0.25 0.25 0.25 0.25 0.25]
  [0.5  0.5  0.5  0.5  0.5 ]
  [0.75 0.75 0.75 0.75 0.75]
  [1.   1.   1.   1.   1.  ]]]


T =
[[[0.         0.         0.         0.         0.        ]
  [0.         0.         0.         0.         0.        ]
  [0.         0.         0.         0.         0.        ]
  [0.         0.         0.         0.         0.        ]
  [0.         0.         0.         0.         0.        ]]

 [[0.78539816 0.78539816 0.78539816 0.78539816 0.78539816]
  [0.78539816 0.78539816 0.78539816 0.78539816 0.78539816]
  [0.78539816 0.78539816 0.78539816 0.78539816 0.78539816]
  [0.78539816 0.78539816 0.78539816 0.78539816 0.78539816]
  [0.78539816 0.78539816 0.78539816 0.78539816 0.78539816]]

 [[1.57079633 1.57079633 1.57079633 1.57079633 1.57079633]
  [1.57079633 1.57079633 1.57079633 1.57079633 1.57079633]
  [1.57079633 1.57079633 1.57079633 1.57079633 1.57079633]
  [1.57079633 1.57079633 1.57079633 1.57079633 1.57079633]
  [1.57079633 1.57079633 1.57079633 1.57079633 1.57079633]]

 [[2.35619449 2.35619449 2.35619449 2.35619449 2.35619449]
  [2.35619449 2.35619449 2.35619449 2.35619449 2.35619449]
  [2.35619449 2.35619449 2.35619449 2.35619449 2.35619449]
  [2.35619449 2.35619449 2.35619449 2.35619449 2.35619449]
  [2.35619449 2.35619449 2.35619449 2.35619449 2.35619449]]

 [[3.14159265 3.14159265 3.14159265 3.14159265 3.14159265]
  [3.14159265 3.14159265 3.14159265 3.14159265 3.14159265]
  [3.14159265 3.14159265 3.14159265 3.14159265 3.14159265]
  [3.14159265 3.14159265 3.14159265 3.14159265 3.14159265]
  [3.14159265 3.14159265 3.14159265 3.14159265 3.14159265]]]


P =
[[[-3.14159265 -1.57079633  0.          1.57079633  3.14159265]
  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]
  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]
  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]
  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]]

 [[-3.14159265 -1.57079633  0.          1.57079633  3.14159265]
  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]
  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]
  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]
  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]]

 [[-3.14159265 -1.57079633  0.          1.57079633  3.14159265]
  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]
  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]
  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]
  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]]

 [[-3.14159265 -1.57079633  0.          1.57079633  3.14159265]
  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]
  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]
  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]
  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]]

 [[-3.14159265 -1.57079633  0.          1.57079633  3.14159265]
  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]
  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]
  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]
  [-3.14159265 -1.57079633  0.          1.57079633  3.14159265]]]

Each of these have a shape of (5, 5, 5), that is \(5^3\) elements in total. We can use them as an input for our function \(f\) to evaluate it for the whole sphere (all possible combinations of the input coordinates).

def f(r, theta, phi):
    """Transform spherical to cartesian coordinates"""

    return (r * np.sin(theta) * np.cos(phi),
            r * np.sin(theta) * np.sin(phi),
            r * np.cos(theta))

np.shape(f(R, T, P))

(3, 5, 5, 5)

Insider tip: Outer operations

We already spoke about the elementwise array-arithmetics before. A really useful feature of NumPy are outer operations, instead of performing the operation elementwise. Before we will look at this, let’s consider the functions np.inner and np.outer on their own and other kinds of products first, which we can connect to linear algebra again.

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

The inner product of two 1D-arrays (vectors) is the scalar product (tensor dot product).

print(np.inner(a, b))
# Ordinary inner product (sum over elementwise products)

8

print(np.dot(a, b))
print(a @ b)
print(np.tensordot(a, b, axes=1))  # tensor dot product

8
8
8

\(\mathbf{a} \cdot \mathbf{b}\)

You can use np.dot also for a general matrix (vector) products.

np.dot(M, a)

array([10, 14])

\(\mathbf{M} \cdot \mathbf{a}\)

np.dot(a, M.T)

array([10, 14])

For vectors in \(\mathbb{R}^3\) you can compute the cross product, too.

np.cross(a, b)

array([-12,   0,   4])

\(\mathbf{a} \times \mathbf{b}\)

The outer product (tensor product) of two 1D arrays (vectors) is a matrix of all possible elementwise products.

print(np.outer(a, b))

print()
print(np.tensordot(a, b, axes=0))

[[ 0  4  0]
 [ 0  8  0]
 [ 0 12  0]]

[[ 0  4  0]
 [ 0  8  0]
 [ 0 12  0]]

\(\mathbf{a} \otimes \mathbf{b}\)

Now, NumPy also allows us to use the outer/inner mechanism not only for products but for any NumPy operation. Using a .outer for example on the np.substract operation will create an array containing the result of the operation performed on all possible pairs of input elements.

A = np.asarray([1, 3, 5, 7])
B = np.asarray([0, 2.0, 4.0, 8.5])

print("A = ", A)
print("B = ", B)

print("\n")
print("np.subtract(A, B) = ")
print(" ", np.subtract(A,B))

print("\n")
print("np.subtract.outer(A, B) = ")
print("    ", np.subtract.outer(A,B))

A =  [1 3 5 7]
B =  [0.  2.  4.  8.5]


np.subtract(A, B) =
  [ 1.   1.   1.  -1.5]


np.subtract.outer(A, B) =
     [[ 1.  -1.  -3.  -7.5]
 [ 3.   1.  -1.  -5.5]
 [ 5.   3.   1.  -3.5]
 [ 7.   5.   3.  -1.5]]

Fix dimensionality related issues with np.newaxis

When we combine multidimensional arrays with each other, we need to take care that their shape can be matched, i.e. broadcasting must be possible. Let’s have a look at these two arrays A and C:

C = np.array([
    [1, 2, 3],
    [4, 5, 6],
    ])

A = np.array([1, 2, 3])

print(C, np.shape(C))
print("\n")
print(A, np.shape(A))

[[1 2 3]
 [4 5 6]] (2, 3)


[1 2 3] (3,)

We can substract A from C without problem.

C - A

array([[0, 0, 0],
       [3, 3, 3]])

Running the cell below will give an error. That is because we want to perform an operation (a subtraction) with C.T and A, having shapes (3, 2) and (3,), respectively.

# C.T = C transposed
C.T - A

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-45-8a6275f83609> in <module>
      1 # C.T = C transposed
----> 2 C.T - A

ValueError: operands could not be broadcast together with shapes (3,2) (3,)

print(C.T, np.shape(C.T))

[[1 4]
 [2 5]
 [3 6]] (3, 2)

We are somehow trying to combine apples and oranges here, but NumPy provides a fix. np.newaxis can be used to extend an array in a dimension. In this case we would like A to be of shape (3, 1), a one dimensional array in two dimensions so to say. Lets try:

C.T - A[np.newaxis]

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-47-eb3ebb4fbf5a> in <module>
----> 1 C.T - A[np.newaxis]

ValueError: operands could not be broadcast together with shapes (3,2) (1,3)

print(A[np.newaxis], np.shape(A[np.newaxis]))

[[1 2 3]] (1, 3)

This failed, too because the new axis (or new dimension) is in the wrong place! We want to conserve the first dimension and append the new one thereafter.

C.T - A[:, np.newaxis]

array([[0, 3],
       [0, 3],
       [0, 3]])

print(A[:, np.newaxis], np.shape(A[:, np.newaxis]))

[[1]
 [2]
 [3]] (3, 1)

Finally, our operation was succesfull.

Note: Since C is a 2D array, we can think of it as a matrix. Using C.T is the transpose of this matrix. In higher dimensions array.T reverses the order of all dimensions.