Matrices and linear algebra in Python

Matrices in Python

In Python, matrices are not it's own thing, but rather a list of list/nested lists. Let's look at an example:

Our matrix of choice will be To represent this matrix in Python, we will consider the matrix as two independent lists,  and , put together in one list. Written in code, it looks like this:

A = [[1, 2, 3], [4, 5, 6]]
# To make the matrices easier to read while coding, we can place the different lists vertical to eachother instead of horizontally.
B = [[1, 2, 3],
	   [4, 5, 6]]


# Printing both of these matrices will result in the same output:
[[1, 2, 3], [4, 5, 6]]

# To alter the output to a more readable format, for-loops are very helpful:
[[1, 2, 3], 
[4, 5, 6]]

Now, let's see how to obtain different values from our matrix:

A = [[1, 2, 3], 
	[4, 5, 6]]

print("A =", A) 				# Entire matrix
print("A[1] =", A[1])      		# 2nd row
print("A[1][2] =", A[1][2])   	# 3rd element of 2nd row
print("A[0][-1] =", A[0][-1])   # Last element of 1st Row

column = [];        			# Empty list
for row in A:
  column.append(row[2])   		# Adding the element in the 3rd column of every row.

print("3rd column =", column)

The script above will result in the following output:

A = [[1, 2, 3], [4, 5, 6]]
A[1] = [4, 5, 6]
A[1][2] = 6
A[0][-1] = 3
3rd column = [3, 6]

Matrices in Python using NumPy

First of all, we have to import the NumPy library. For the sake of this tutorial, we will only do it once, but you can assume it is done everywhere NumPy is being used.

import numpy as np

To make an array (or matrix) using NumPy, we will use the function numpy.array, and simply use the same syntax as before, but now as a function parameter.

A = np.array([[1, 2, 3], [4, 5, 6]])		# Array of ints

B = np.array([[1.1, 2, 3], [4, 5, 6]])		# Array of floats

C = np.array([[1, 2, 3], [4, 5, 6]], dtype=complex)		# Array of complex numbers

When printing the matrices above using print, the output will look like this:

[[1 2 3]					# A
 [4 5 6]]

[[1.1 2.  3. ]				# B
 [4.  5.  6. ]]

[[1.+0.j 2.+0.j 3.+0.j]		# C
 [4.+0.j 5.+0.j 6.+0.j]]

Notice how the matrices are printed in an easily readable format without the use of for-loops. This is another benefit of using NumPy arrays.

There are several other functions in NumPy that can create specific matrices. Here are some of them:

• numpy.zeros returns a new array with given shape and type, filled with zeros.

  zeros Expand source

  A = np.zeros((2, 3))
  
  # When printed:
  
  [[0. 0. 0.]
   [0. 0. 0.]]

• numpy.ones returns a new array with given shape and type, filled with ones.

  ones Expand source

  A = np.ones((2, 3))
  
  # When printed:
  
  [[1. 1. 1.]
   [1. 1. 1.]]

• numpy.empty returns a new array with given shape and type, without initalizing entries.

  empty Expand source

  A = np.empty((2, 3))
  
  # When printed:
  
  [[ 1.33360294e+241  1.00200641e-310  3.42196482e-210]		# Random numbers
   [ 2.42869050e-313  1.59166660e-308 -1.95565444e-310]]

• numpy.full returns a new array with given shape and type, filled with fill_value.

  full Expand source

  A = np.full((2, 3), 2019)
  
  B = np.full((2, 3), np.inf)
  
  # When printed:
  
  [[2019 2019 2019]		# A
   [2019 2019 2019]]	
  
  [[inf inf inf]			# B
   [inf inf inf]]

Linear algebra using NumPy

NumPy has several useful functions for built-in linear algebra. For a full list, check out the NumPy documentation on linear algebra. Optionally, and as mentioned at the top of this page, see SciPys documentation on linear algebra.

Dot product/Matrix multiplication

To compute the dot product of two arrays, A and B, we can use the function numpy.dot. However, as recommended by NumPy, if both arrays are 2-D arrays it is preffered to use numpy.matmul or "A @ B". As an example, we will now compute the dot product of two matrices, A and B, in all three ways. Let 

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

# Using "np.dot"
print(np.dot(A, B))

# Using "np.matmul"
print(np.matmul(A, B))

# Using "@"
print(A @ B)


# All three version will produce the same output:
[[19 22]
 [43 50]]

Solving a system of linear equations

Let the set of equations we would like to solve be and . Let's first write it as matrices:

The function we will use to solve the equation is numpy.linalg.solve. It takes two parameters, a and b, where both are arrays (remember that our matrix is an array type). numpy.linalg.solve will then computes the “exact” solution, x, of the well-determined, i.e., full rank, linear matrix equation ax = b .

a = np.array([[5, 2],[1,4]])    # Left-hand side
b = np.array([33, 21])          # Right-hand side

print(np.linalg.solve(a,b))     # Printing the solution


# Output:
[5. 4.]		# In our example this evaluates to x = 5, y = 4, which is the correct answer!


# You may also save the values directly into variables:
x, y = np.linalg.solve(a,b)

Eigenvalues and eigenvectors

numpy.linalg.ein is a function that computes the eigenvalues and eigenvectors of a square array, returning two arrays consisting of the eigenvalues and eigenvectors respectivly. Lets look at an example. We want to find the eigenvalues and eigenvectors of matrix A, where.

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

w, v = np.linalg.eig(A)		# "w" are the eigenvalues, "v" the eigenvectors

print(w)
print(v)

 
# Output:
[ 3. -5.]						# w
[[ 0.83205029 -0.4472136 ]		# v
 [ 0.5547002   0.89442719]]

The eigenvectors are arranged so that column v[:,i] is the eigenvector corresponding to the eigenvalue w[i], i.e. is the eigenvector to the eigenvalue 3, and is the eigenvector to the eigenvalue -5.

Note that the function returns the eigenvectors with a unit length. An easier way to present the eigenvectors above is as and , respectively.

Determinant

To compute the determinant of a square array, we will use the function numpy.linalg.det. It takes an array (or several arrays combined) as input, and returns the determinant. As an example, we will compute the determinant of matrix A, where.

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

print(np.linalg.det(A))


# Output:
-2.0000000000000004

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

 # To compute the determinant of multiple arrays at once, we have to put our arrays into a single array as input.
 # This will result in an array of length 3 as output, with the determinants respectivly.
print(np.linalg.det([A, B, C]))    


# Output:
[-2. -3.  4.] 		# Evaluates to: det(A) = -2, det(B) = -3 and det(C) = 4.

Inverting a matrix

To invert a matrix, we can use the function numpy.linalg.inv. It takes a square array as input, and returns the invers of the given array. As an example, we will compute the invers of matrix A, where 

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

print(np.linalg.inv(A))


# Output:
[[-2.   1. ]
 [ 1.5 -0.5]]

# To check if the invers we got in the example above is correct, we simply compute the dot product of A and its invers.
A = np.array([[1, 2], [3, 4]])
B = np.array([[-2., 1. ], [1.5, -0.5]])		# This is the presumed invers of A.

print(A @ B)	# Matrix product/ dot product


# Output:
[[1. 0.]		# This is the identity matrix, which implies that B really is the invers of A.
 [0. 1.]]		# Remember how the invers works, where "B @ A" would have given us the same answer.

Exercises and solutions

If you want have a look at some exercises and solutions regarding NumPy, SciPy, matrices and linear algebra, visit Exercises and solutions, linear algebra in Python.