Linear Algebra with NumPy

Installing NumPy¶

In many mathematical applications, calculations with vectors and matrices are required, for example in the numerical computation of integrals, differential equations, or problems from graph theory.

In this chapter, we will focus on a Python library that provides appropriate data types for vectors and matrices, as well as efficient operations on them: NumPy.

First, the NumPy library must be installed. When using Conda, this is done with the console command

conda install numpy

Next, the library must be imported into our Python script. This could – as in section Variables and Data Types – be done with

from numpy import *

However, this approach is not recommended, since NumPy also provides functions like sin, cos, sqrt, etc., which would overwrite identically named functions from other libraries (e.g., math).

Instead, we typically use the following notation:

import numpy as np

The addition as np simply indicates that we can refer to the library in the future by the shorter name np instead of the longer numpy.

Working with Vectors¶

Creating Vectors¶

A vector or NumPy array can be initialized using the function np.array(...) from any iterable object (list, tuple, ...), whose elements are of the same type:

a = np.array([1.,2.,3.])
b = np.array((6.,5.,4.))
print("a is a", type(a), "with value", a)
print("b is a", type(b), "with value", b)

a is a <class 'numpy.ndarray'> with value [1. 2. 3.]
b is a <class 'numpy.ndarray'> with value [6. 5. 4.]

The class ndarray represents a multidimensional array. In our case, a vector.

The number of dimensions can be obtained using

a.ndim

1

The data type of the stored elements is

a.dtype

dtype('float64')

The shape of the array (i.e., the size of each dimension) can be obtained using

a.shape

(3,)

Two other useful functions for creating NumPy arrays are linspace and arange:

x = np.linspace(0,3,6) # Equally spaced grid for [0,3] with 6 points
x

array([0. , 0.6, 1.2, 1.8, 2.4, 3. ])

The function np.linspace(start, stop, num) generates num evenly spaced points in the interval [start, stop]. By default, both endpoints are included.

x = np.arange(0, 3, 0.5) # Equally spaced grid for [0,3) with increment 0.5
x

array([0. , 0.5, 1. , 1.5, 2. , 2.5])

The function np.arange(start, stop, step) generates a sequence of values starting at start with a step size of step. The end value stop is not reached (similar to the Python range function).

Basic Vector Operations¶

An ndarray is – similar to a list or a tuple – an iterable object. Therefore, we can iterate over all elements using a for loop:

val = 0
for e in a:
    val += e
print("The sum of the entries in", a, "is", val)

The sum of the entries in [1. 2. 3.] is 6.0

With NumPy arrays, basic arithmetic operations can also be performed directly:

a+b

array([7., 7., 7.])

b-a

array([5., 3., 1.])

a*b

array([ 6., 10., 12.])

a/b

array([0.16666667, 0.4 , 0.75 ])

We observe that these operations are performed element-wise.

For addition and subtraction, this exactly corresponds to the definition of vector addition from linear algebra. For multiplication and division, this is less obvious – one might have expected the dot product or the cross product instead. We will learn about these operations later.

The arithmetic operations +, -, *, and / only work if the involved arrays have the same shape. Otherwise, an error occurs:

c = np.array([8,7,6,5])
a+c

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[13], line 2
      1 c = np.array([8,7,6,5])
----> 2 a+c

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

Now let’s examine which methods the ndarray class provides. To do this, we can type a.<TAB> in a notebook or call the function dir(a). In both cases, we get a long list of available methods.

Let’s test some of these methods:

print("The smallest entry of a is", a.min(), "at index", a.argmin())
print("The dot product <a,b> is", a.dot(b))
print("The sum of all entries of a is", a.sum())

The smallest entry of a is 1.0 at index 0
The dot product <a,b> is 28.0
The sum of all entries of a is 6.0

Additional Arithmetic Operations¶

In addition to the basic methods directly provided by the ndarray class, there are many other arithmetic operations implemented as free functions in the numpy library.

If we type np.<TAB>, we get a list of these functions. Notice that functions like sqrt, exp, sin, or cos are also defined here. Although these exist in the math library, they cannot be applied to NumPy arrays there:

import math
math.exp(a)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[15], line 2
      1 import math
----> 2 math.exp(a)

TypeError: only 0-dimensional arrays can be converted to Python scalars

This results in an error because math.exp is only defined for scalar values.

The corresponding function from the numpy library, however, is applied element-wise to the array:

np.exp(a)

array([ 2.71828183, 7.3890561 , 20.08553692])

The exponential function is applied to each component of the vector.

Other functions such as log, sin, cos, tan, abs, etc., are also applied element-wise to NumPy arrays.

Similarly, comparison operations work element-wise:

x = np.array([0,1,2,3,4])
y = (x <=2)
y

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

The result is an array of boolean values, indicating for each element of x whether the condition is satisfied.

In addition to these element-wise operations, the numpy module also provides vector-specific operations, such as the dot product or the cross product.

np.dot(a,b)

np.float64(28.0)

np.cross(a,b)

array([-7., 14., -7.])

Upon closer inspection of the functions in the numpy module, one finds that many operations from linear algebra are organized in the submodule numpy.linalg. These include, for example, functions for computing norms, determinants, or matrix inverses.

If we type

np.linalg.<TAB>

we get a list of available functions.

The function np.linalg.norm can be used to compute various vector norms:

print("Euclidean norm :", np.linalg.norm(a))
print("Maximum norm   :", np.linalg.norm(a, np.inf))
print("1-norm         :", np.linalg.norm(a, 1))

Euclidean norm : 3.7416573867739413
Maximum norm   : 3.0
1-norm         : 6.0

Accessing Vector Elements¶

Finally, let’s look at how to access individual or multiple elements of a NumPy array. This is done – just like with lists – using the [] operator:

a = np.linspace(0,1,11)
a

array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ])

print("Fourth entry:") # Index counting starts at 0
a[3]

Fourth entry:

np.float64(0.30000000000000004)

We can also extract a part of the array by specifying an index range i:j:

a[3:6] # indices from 3 to 5

array([0.3, 0.4, 0.5])

In general:

i is the start index (inclusive)
j is the end index (exclusive)

More examples:

a[:6] # indices up to 5

array([0. , 0.1, 0.2, 0.3, 0.4, 0.5])

a[6:] # indices from 6 onward

array([0.6, 0.7, 0.8, 0.9, 1. ])

We can also overwrite sub-vectors, of course respecting the dimensions:

a[1:4] = np.ones((3,))  # Write ones into entries 1 to 3
a[6:8] = np.zeros((2,)) # Write zeros into entries 6 and 7
a

array([0. , 1. , 1. , 1. , 0.4, 0.5, 0. , 0. , 0.8, 0.9, 1. ])

In the process, we have also learned two more functions that can create arrays with constant entries:

np.ones(shape) creates an array of the given shape, where all entries are 1.
np.zeros(shape) creates an array of the given shape, where all entries are 0.

Working with Matrices¶

Creating Matrices¶

For matrices, we also use the data type numpy.ndarray. To initialize, we pass a list of lists to the function np.array, where each inner list corresponds to a row of the matrix. NumPy automatically recognizes that a two-dimensional array, i.e., a matrix, should be created:

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

print("A is of type:", type(A))
print("A has shape:", A.shape)
print(A)

A is of type: <class 'numpy.ndarray'>
A has shape: (3, 3)
[[1. 4. 2.]
 [2. 0. 3.]
 [1. 1. 4.]]

The shape of the matrix (A.shape) gives the number of rows and columns, i.e., (3,3) in this example.
Internally, A is still an ndarray, just with two dimensions (A.ndim == 2).
Just like with vectors, we can access individual elements or subarrays, e.g., A[0,1] for row 0, column 1.

There are several ways to create matrices in NumPy:

Identity Matrix
A square identity matrix can be created directly using np.eye(n):

B = np.eye(3)
B

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

Empty Matrix and Manual Filling
With np.zeros((m,n)) a matrix of size m×n can be created, with all entries initially set to 0. Individual entries can then be assigned directly:

C = np.zeros((3,3)) # The argument is a tuple specifying the size of the matrix
C[0,0] = 1
C[1,1] = 2
C[1,2] = 5
C[2,1] = 4
C[2,2] = 1
C

array([[1., 0., 0.],
       [0., 2., 5.],
       [0., 4., 1.]])

Diagonal Matrix from a Vector
From a given vector v, a diagonal matrix can be created where the elements of v appear on the main diagonal:

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

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

Matrix Arithmetic¶

For matrices, the basic arithmetic operations work similarly to vectors: addition, subtraction, multiplication, and division, as well as many functions from the numpy package (exp, sin, cos, log, etc.) are applied element-wise to the matrices.

C = A+B
C

array([[2., 4., 2.],
       [2., 1., 3.],
       [1., 1., 5.]])

D = A*B
D

array([[1., 0., 0.],
       [0., 0., 0.],
       [0., 0., 4.]])

E = np.exp(A)
E

array([[ 2.71828183, 54.59815003,  7.3890561 ],
       [ 7.3890561 ,  1.        , 20.08553692],
       [ 2.71828183,  2.71828183, 54.59815003]])

Additional operations are provided by the numpy.ndarray class, such as transposing a matrix:

C.transpose()

array([[2., 2., 1.],
       [4., 1., 1.],
       [2., 3., 5.]])

For linear algebra operations, one uses the functions from numpy or numpy.linalg:

Matrix Multiplication:

np.matmul(A,B) # Free function
A@B            # Alternatively using the @ operator

array([[1., 4., 2.],
       [2., 0., 3.],
       [1., 1., 4.]])

Determinant:

np.linalg.det(C)

np.float64(-22.000000000000004)

Inverse:

np.linalg.inv(C)

array([[-0.09090909,  0.81818182, -0.45454545],
       [ 0.31818182, -0.36363636,  0.09090909],
       [-0.04545455, -0.09090909,  0.27272727]])

Eigenvalues- and vectors:

E,V = np.linalg.eig(C)
print("Eigenvalues    :\n", E)
print("Eigenvectors   :\n", V)

Eigenvalues    :
 [ 6.97413644 -1.33574651  2.36161007]
Eigenvectors   :
 [[ 0.63932126  0.77319109 -0.8514755 ]
 [ 0.50515119 -0.63379132 -0.29406667]
 [ 0.57973321 -0.02200211  0.43418229]]

Note: The eigenvectors are stored column-wise in V. Let’s test the calculation by checking $C\,v - \lambda\,v=0$ :

for i in range(3):
    res = np.matmul(C, V[:,i]) - E[i]*V[:,i]  # C*v-lambda*v    
    print("Error: ", np.linalg.norm(res))

Error:  3.66205343881779e-15
Error:  1.5852031831917945e-15
Error:  2.434909185329438e-15

Solving Linear Systems of Equations:

b = np.array([24.,23.,30.])
x = np.linalg.solve(C, b)
print("x =", x)

# Test:
np.matmul(C,x)-b

x = [3. 2. 5.]

array([ 0.00000000e+00, 0.00000000e+00, -3.55271368e-15])

Accessing Matrix Entries

Accessing individual entries or submatrices works similarly to one-dimensional ndarray arrays:

print("Entry C_11           :", C[0,0])
print("Second column        :", C[:,1])
print("Third row            :", C[2,:])
print("Last (=third) row    :", C[-1,:])

Entry C_11           : 2.0
Second column        : [4. 1. 1.]
Third row            : [1. 1. 5.]
Last (=third) row    : [1. 1. 5.]

Stacking Matrices

With the functions numpy.vstack and numpy.hstack, matrices can be concatenated vertically and horizontally, respectively:

np.vstack([A,B,C]) # Stacks A, B, C vertically

array([[1., 4., 2.],
       [2., 0., 3.],
       [1., 1., 4.],
       [1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.],
       [2., 4., 2.],
       [2., 1., 3.],
       [1., 1., 5.]])

np.hstack([A,B,C]) # Stacks A, B, C horizontally

array([[1., 4., 2., 1., 0., 0., 2., 4., 2.],
       [2., 0., 3., 0., 1., 0., 2., 1., 3.],
       [1., 1., 4., 0., 0., 1., 1., 1., 5.]])

If you want to “attach” a vector to a matrix, the vector may first need to be reshaped into an $n \times 1$ matrix:

np.vstack([A,b])

array([[ 1.,  4.,  2.],
       [ 2.,  0.,  3.],
       [ 1.,  1.,  4.],
       [24., 23., 30.]])

np.hstack([A,b.reshape((3,1))])

array([[ 1.,  4.,  2., 24.],
       [ 2.,  0.,  3., 23.],
       [ 1.,  1.,  4., 30.]])

Note: reshape((n,1)) reshapes the vector into a column matrix so that hstack works correctly.

Mutable or Immutable?

NumPy arrays behave as mutable objects. This means that changes within functions or through new variables can modify the original array:

def modify_matrix(A):
    A[1,1] = 1.
    
A = np.zeros((3,3))
modify_matrix(A)
print(A)

[[0. 0. 0.]
 [0. 1. 0.]
 [0. 0. 0.]]

The array A was modified by the function since ndarray is mutable.

Let us consider a simple assignment:

B = A
B[1,1] = 2.
A

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

Here, the name B refers to the same object as A. Therefore, changes made via B directly affect A.

If you want to create a true copy that is independent of the original, use:

B = np.copy(A)
B[1,1] = 3.
print(A)

[[0. 0. 0.]
 [0. 2. 0.]
 [0. 0. 0.]]

A change to B does not affect A, since np.copy creates a new array.

NumPy Glossary – Vectors & Matrices¶

Creating Arrays¶

Function	Description	Example
`np.array(seq)`	Converts an iterable object into an array	`np.array([1,2,3])`
`np.zeros(shape)`	Array filled with zeros	`np.zeros((3,3))`
`np.ones(shape)`	Array filled with ones	`np.ones(5)`
`np.eye(n)`	Identity matrix (n×n)	`np.eye(3)`
`np.diag(v)`	Diagonal matrix from vector v	`np.diag([1,2,3])`
`np.linspace(start, stop, num)`	Evenly spaced values	`np.linspace(0,1,5)`
`np.arange(start, stop, step)`	Values with fixed step size	`np.arange(0,3,0.5)`

Basic Attributes of Arrays (`ndarray`)¶

Attribute	Meaning
`.ndim`	Number of dimensions
`.shape`	Shape / size of the array
`.size`	Number of elements
`.dtype`	Data type of the elements

Element-wise Operations¶

Operator / Function	Description
`+,-,*,/`	Element-wise addition, subtraction, multiplication, division
`np.exp(a)`	Element-wise exponential function
`np.log(a)`	Element-wise natural logarithm
`np.sin(a), np.cos(a), np.tan(a)`	Element-wise trigonometric functions
Comparison operators (`<`, `<=`, `>`, `>=`, `==`, `!=`)	Element-wise

Methods of `ndarray`¶

Method	Description
`.sum()`	Sum of all elements
`.min()`	Smallest value
`.max()`	Largest value
`.argmin()`	Index of the smallest value
`.argmax()`	Index of the largest value
`.dot(b)`	Dot product (vectors only)
`.transpose()`	Transposed matrix

Free Functions in `numpy`¶

Function	Description
`np.dot(a,b)`	Dot product / matrix-vector product
`np.cross(a,b)`	Cross product (3D vectors)
`np.matmul(a,b)`	Matrix multiplication
`np.vstack([A,B])`	Vertical stacking of matrices
`np.hstack([A,B])`	Horizontal stacking of matrices

Functions in `numpy.linalg` (Linear Algebra)¶

Function	Description
`np.linalg.norm(a, ord=2)`	Vector norm (default: 2-norm)
`np.linalg.inv(A)`	Inverse of a matrix
`np.linalg.det(A)`	Determinant
`np.linalg.eig(A)`	Eigenvalues and eigenvectors
`np.linalg.solve(A,b)`	Solves linear system Ax=b

Accessing Elements¶

Syntax	Meaning
`a[i]`	i-th entry of a vector
`a[i:j]`	entries from i to j-1
`a[:j]`	entries from 0 to j-1
`a[i:]`	entries from i to end
`A[i,j]`	entry (i,j) of a matrix
`A[i,:]`	row i
`A[:,j]`	column j
`A[i:j, k:l]`	submatrix

Mutable Behavior¶

NumPy arrays are mutable. Changes within functions or through assignments affect the original array.
For true copies: B = np.copy(A)

Installing NumPy¶

Working with Vectors¶

Creating Vectors¶

Basic Vector Operations¶

Additional Arithmetic Operations¶

Accessing Vector Elements¶

Working with Matrices¶

Creating Matrices¶

Matrix Arithmetic¶

NumPy Glossary – Vectors & Matrices¶

Creating Arrays¶

Basic Attributes of Arrays (ndarray)¶

Element-wise Operations¶

Methods of ndarray¶

Free Functions in numpy¶

Functions in numpy.linalg (Linear Algebra)¶

Accessing Elements¶

Mutable Behavior¶

Basic Attributes of Arrays (`ndarray`)¶

Methods of `ndarray`¶

Free Functions in `numpy`¶

Functions in `numpy.linalg` (Linear Algebra)¶