Numpy Array Operations and Linear Algebra¶

Note: We will see more tools for linear algebra in the section on Scipy Tools, which introduces the package Scipy.

As of version 3.5, Python can handle matrix multiplication on Numpy arrays, using the at sign “@”:

C = A @ B

Why this strange notation? Because

D = A * B

is instead “point-wise” multiplication of arrays, with D[i,j] = A[i,j] * B[i,j] instead of

\[C_{i,j}= \sum_{k} a_{ik} b_{kj}\]

Aside 1: Numpy matrix is now redundant!

Previously, Numpy supported matrix calculations withe variables of type matrix, a special sub-tpye of array. However, with the addition of “@” for matrix multiplication of arrays, that option is now un-needed, and I advise you to avoid it.

Aside 2: Floating point numbers vs integers

Though Python is generally good at understanding when an integer like 7 is to be used as a floating point (real) number, it is sometimes best to make this distinction explicitly when working with module numpy; otherwise sometimes division done within numpy functions returns an integer answer, like 7/2 = 3.

Thus from now on, when I mean an integer to be used as a floating point number, I give it a decimal point: 7./2. will reliably be 3.5

First we explore some basic features offered by Numpy; then we will look at some more advanced tools from package Scipy, and specifically its module scipy.linalg for linear algebra.

import numpy as np

A = np.array([[1, 2, 3],[4, 5, 6]])
B = np.array([[1, 2], [3,4], [5, 6]])
C = np.array([[10, 9, 8],[7, 6, 5]])

print(f"A is\n{A}")
print(f"B is\n{B}")
print(f"C is\n{C}")

A is
[[1 2 3]
 [4 5 6]]
B is
[[1 2]
 [3 4]
 [5 6]]
C is
[[10  9  8]
 [ 7  6  5]]

print(f"The matrix product A times B is:\n {A @ B}")
print(f"The matrix product B times A is:\n {B @ A}")

The matrix product A times B is:
 [[22 28]
 [49 64]]
The matrix product B times A is:
 [[ 9 12 15]
 [19 26 33]
 [29 40 51]]

print(f"The array product A * B fails:\n{A * B}")

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-4-3bb26336abf9> in <module>
----> 1 print(f"The array product A * B fails:\n{A * B}")

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

On the other hand:

print(f"Array product A times C is\n{A * C}\n")

Array product A times C is
[[10 18 24]
 [28 30 30]]

print(f"Matrix product A times C fails:\n{A @ C}")

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-6-55b834c682c2> in <module>
----> 1 print(f"Matrix product A times C fails:\n{A @ C}")

ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 2 is different from 3)

The matrix transpose \(A^T\) is given by `A.T`¶

print(f"A-transpose is\n{A.T}")

A-transpose is
[[1 4]
 [2 5]
 [3 6]]

Slicing: Extracting rows, columns, and other rectangular chunks from matrices¶

This works with Python lists and Numpy arrays, and we have seen some of it before; I review it here because it will help with doing the row operations of linear algebra.

Index notation for slicing¶

For an index with n possible values, from 0 to n-1:

a:b means indices i in th usual semi-open interval \(a \leq i < b\) or \([a,b)\)
a: is short for a:n, so indices \(a \leq i\), all the way to the maximum index value
:b is short for 0:b, so all indices \(i < b\)
: combines both of the above, so is short for 0:n; all possible indices
index value -1 refers the last entry; the same as index n-1, but without needing to know n.
index value -m refers the “m-th last” entry; the same as index n-m, again without needing to know n.

print(f"A:\n{A}")
print(f"The column of index 1 (presented as a row vector): {A[:,1]}")
print(f"The row of index 1: {A[1,:]}")
print(f"The first 2 elements of the row of index 1: {A[1,:2]}")
print(f"Another way to say the above: {A[1,0:2]}")
print(f"The bottom-right element: {A[-1,-1]}")
print(f"The 2x2 sub-matrix in the bottom-right corner:\n{A[-2:,-2:]}")

A:
[[1 2 3]
 [4 5 6]]
The column of index 1 (presented as a row vector): [2 5]
The row of index 1: [4 5 6]
The first 2 elements of the row of index 1: [4 5]
Another way to say the above: [4 5]
The bottom-right element: 6
The 2x2 sub-matrix in the bottom-right corner:
[[2 3]
 [5 6]]

Synonyms for array names with the equal sign (not copying!)¶

If we use the equal sign between two array names, that makes them synonyms, referring to the same values:

A = np.array([[1, 2, 3],[4, 5, 6]])
print(f"A is\n{A}")
Anickname = A
print(f"Anickname is\n{Anickname}")
Anickname[0,0] = 12
print(f"Anickname is now\n{Anickname}")
print(f"and so is A!:\n{A}")

A is
[[1 2 3]
 [4 5 6]]
Anickname is
[[1 2 3]
 [4 5 6]]
Anickname is now
[[12  2  3]
 [ 4  5  6]]
and so is A!:
[[12  2  3]
 [ 4  5  6]]

Copying arrays with method .copy()¶

Thus if we want a separate new array or list with the same elements initially, we must make a copy with the method .copy(), not the equal sign:

A = np.array([[1, 2, 3],[4, 5, 6]])
print(f"A is\n{A}")
Acopy = A.copy()
print(f"Acopy is\n{Acopy}")
Acopy[0,0] = 54
print(f"Acopy is now\n{Acopy}")
print(f"A is still\n{A}")

A is
[[1 2 3]
 [4 5 6]]
Acopy is
[[1 2 3]
 [4 5 6]]
Acopy is now
[[54  2  3]
 [ 4  5  6]]
A is still
[[1 2 3]
 [4 5 6]]

Exercise A¶

Create a Numpy array (not a Numpy matrix; those are now mostly obsolete!) containing the matrix

\[\begin{split}A = \left[ \begin{array}{ccc} 4. & 2. & 1. \\ 9. & 3. & 1. \\ 25. & 5. & 1. \end{array} \right]\end{split}\]

and one containing the vector

\[b = [ 0.693147, 1.098612, 1.386294]\]

Exercise B¶

Create arrays \(c\) and \(d\) containing respectively the last row of \(A\) and the middle column of \(A\).

Note: Do this by manipulating the array A with indexing and slicing operations, without entering any numerical values for array entries explicity.

Python for Scientific Computing