Introduction To Matrices| Linear Algebra Using Python

Linear algebra consists of the methods to solve the system of linear equations and their applications to tackle real-world problems. This system of equations is expressed in a regular fashion to create an array or a ‘matrix’ of elements. These matrices allow us to perform different operations on them and make it easy to solve the system of equations.

In this article, we will understand what is a matrix and how to represent a system of linear equations in the form of a matrix. We will also go through various matrix creation routines in Python using the NumPy library.

What is a Matrix
Representing System of Linear Equations with Matrices
Representation of Matrix Elements
Creating Matrices in Python using NumPy

What is a Matrix

A matrix is an ordered rectangular array of numbers (real or complex), which has m rows (horizontal set of elements) and n columns (vertical set of elements). In general, rows and columns are called lines. The order of a matrix is represented as m x n (pronounced m by n).

Elements of the matrix are enclosed in [ ] or ( ). For example-

$$\begin{aligned}
A =
\begin{bmatrix}
5 & 3 & 58 \\
-4 & 23 & 11 \\
34 & 2 & -67
\end{bmatrix}
\end{aligned}$$

Representing System of Linear Equations with Matrices

Let’s now see how we can represent a system of linear equations in the form of a matrix. The equations contain the variable with coefficients and a constant term.

$$\begin{aligned}
2x + 5y &= 16\\
4x + 3y &= 8
\end{aligned}$$

These equations are represented with an augmented matrix, where the first equation forms the first row and the second equation forms the second row.

The first column represents the $x$ coefficients and the second column represents the $y$ coefficients. The constants in the equations are entered in the third column.

The augmented matrix formed from the above system of equations, say A, is:

$$\begin{aligned}
A =
\begin{bmatrix}
2 & 5 & 16 \\
4 & 3 & 8
\end{bmatrix}
\end{aligned}$$

If a term is missing in the equation, it means that the coefficient of that variable is 0.

$$\begin{aligned}
8y – 3z &= 19\\
2x + z &= -5\\
7x + 4y + 6z &= 28\\
\Downarrow\\
0x + 8y – 3z &= 19\\
2x + 0y + z &= -5\\
7x + 4y + 6z &= 28\\
\end{aligned}$$ $$\begin{aligned}
\Rightarrow B =
\begin{bmatrix}
0 & 8 & -3 & 19 \\
2 & 0 & 1 & -5 \\
7 & 4 & 6 & 28
\end{bmatrix}
\end{aligned}$$

If the constant terms are on the left side of the equal sign or variable terms are on the right side of the equal sign, we have to rearrange the equations to bring them to the standard form before representing them in the augmented matrix.

Also, make sure that the sequence of the variable terms is uniform.

$$\begin{aligned}
3 – 5x + 7z &= 12y\\
4x + y – 11 &= -9z\\
6y + 8x &= 0\\
\Downarrow\\
– 5x -12y + 7z &= -3\\
4x + y + 9z &= 11\\
8x + 6y + 0z &= 0\\
\end{aligned}$$ $$\begin{aligned}
\Rightarrow C =
\begin{bmatrix}
-5 & -12 & 7 & -3 \\
4 & 1 & 9 & 11 \\
8 & 6 & 0 & 0
\end{bmatrix}
\end{aligned}$$

Representation of Matrix Elements

Consider the matrix $A$ with 2 rows and 4 columns. An individual element of the matrix is represented as $a_{ij}$, where $i$ is the row number and $j$ is the column number. The entire 2 x 4 matrix can be represented as:

$$\begin{aligned}
A_{2\times4} =
\begin{bmatrix}
a_{11} & a_{12} & a_{13} & a_{14} \\
a_{21} & a_{22} & a_{23} & a_{24}
\end{bmatrix}
\end{aligned}$$

Creating Matrices in Python using NumPy

Now we will see the ways to create the matrices in Python using the NumPy library.

If you don’t have the NumPy library installed on your computer, use the Python package manager (pip) to get it installed. Type the following line in your command prompt (terminal):

pip install numpy

Check my article on setting up Python for science to get started.

To use the functions available in the NumPy library, we first need to import it using:

import numpy as np

Here np is called an alias to numpy and is used to avoid typing the word numpy every time we use any function from the library.

One thing to mention before we dive in is that the matrices can be made up of integers, floating-point numbers, complex numbers, or boolean values. These are called the data types in the programming world.

Creating the matrices from the existing data

The simplest way to create a matrix from the existing data is to use the numpy.array( ) function. Let’s consider that we have the data in the form of a list or a tuple. The numpy.array( ) function takes mainly two arguments, first is the list (or tuple) and second is the data type (dtype). The second argument is optional and the data type will be detected from the list provided by us if not specified. For example:

A = np.array([1,2,3])
print('Row Matrix of integers:')
print(A)

B = np.array([[1.0,2.0,3.0],[4.0,5.0,6.0]])
print('Rectangular Matrix of floats:')
print(B)

Output:

Row Matrix of integers:
[1 2 3]
Rectangular Matrix of floats:
[[1. 2. 3.]
 [4. 5. 6.]]

Here the row matrix and the rectangular matrix are two of the types of matrices.

The data type will be enforced if we specify it using the dtype argument to the array function.

my_list = [1,2,3,4,5]

P = np.array(my_list,dtype=int)
print('Matrix P with integer elements:')
print(P)

my_tuple = ((1,2,3),(4,5,6))

Q = np.array(my_tuple ,dtype=float)
print('Matrix Q with real number elements:')
print(Q)

R = np.array(my_tuple ,dtype=complex)
print('Matrix R with complex number elements:')
print(R)

Output:

Matrix P with integer elements:
[1 2 3 4 5]
Matrix Q with real number elements:
[[1. 2. 3.]
 [4. 5. 6.]]
Matrix R with complex number elements:
[[1.+0.j 2.+0.j 3.+0.j]
 [4.+0.j 5.+0.j 6.+0.j]]

Another function we can use to create a Row Matrix from an existing iterable is the numpy.fromiter( ) function.

This function takes three arguments. First is the iterable, second is the dtype and count is an optional argument. count is used to specify the number of items to be read from the iterable with the default value of -1, i.e., read all the elements.

my_list = [1,2,3,4,5] 
iter_1 = iter(my_list)
iter_2 = (x*x*x for x in range(1,10))

A = np.fromiter(iter_1,float)
B = np.fromiter(iter_2,int,5)

print('Matrix A:')
print(A)
print('Matrix B with elements as cubes of first 5 integers:')
print(B)

Output:

Matrix A:
[1. 2. 3. 4. 5.]
Matrix B with elements as cubes of first 5 integers:
[  1   8  27  64 125]

We can also reshape an existing matrix to create a new matrix with a different shape, provided that the total number of elements remains the same. We will use the reshape( ) method associated with the NumPy array.

A = np.array(range(1,13)) # 12 elements

B = A.reshape((3,4))
C = B.reshape((2,6))

print('Row Matrix A:')
print(A)
print('Matrix B:')
print(B)
print('Matrix C:')
print(C)

Output:

Row Matrix A:
[ 1  2  3  4  5  6  7  8  9 10 11 12]
Matrix B(3x4):
[[ 1  2  3  4]
 [ 5  6  7  8]
 [ 9 10 11 12]]
Matrix C(2x6):
[[ 1  2  3  4  5  6]
 [ 7  8  9 10 11 12]]

Along with the in-built range( ) function used above, we can make use of the sequence generation routines from the NumPy library, which include the arange( ), linspace( ) and logspace( ) functions.

The numpy.arange( ) function can be used in the same way as the range( ) function. The only difference is that we can have the real numbers (float data type) as the start, stop and step arguments of the function.

P = np.arange(2.5,4.8,0.2)
Q = P.reshape((3,4))

print('Row Matrix P:')
print(P)
print('Matrix P:')
print(Q)

Output:

Row Matrix P:
[2.5 2.7 2.9 3.1 3.3 3.5 3.7 3.9 4.1 4.3 4.5 4.7]
Matrix P:
[[2.5 2.7 2.9 3.1]
 [3.3 3.5 3.7 3.9]
 [4.1 4.3 4.5 4.7]]

The numpy.linspace( ) and numpy.logspace( ) functions require the start, stop and the number of elements (num) as the primary arguments. In case of the logspace( ) function, we can also specify the base of the log space. The default value of base is 10.

L = np.reshape(np.linspace(1,4,num=9),(3,3))

M = np.reshape(np.logspace(1,4,num=9),(3,3))
N = np.reshape(np.logspace(1,10,num=10,base=2),(2,5))

print('Matrix L:')
print(L)
print('------------------------------')
print('Matrix M:')
print(M)
print('------------------------------')
print('Matrix N:')
print(N)

Output:

Matrix L:
[[1.    1.375 1.75 ]
 [2.125 2.5   2.875]
 [3.25  3.625 4.   ]]
------------------------------
Matrix M:
[[   10.            23.71373706    56.23413252]
 [  133.35214322   316.22776602   749.89420933]
 [ 1778.27941004  4216.96503429 10000.        ]]
------------------------------
Matrix N:
[[   2.    4.    8.   16.   32.]
 [  64.  128.  256.  512. 1024.]]

Creating the matrices with zeros and ones

Various functions are available in NumPy to create matrices with different combinations of zeros and ones.

Null or Zero Matrix

The numpy.zeros( ) function is used to create a Null or Zero Matrix. This function takes the tuple of the shape of a matrix as an input and the data type is an optional argument.

L = np.zeros(3)
M = np.zeros((2,4))
N = np.zeros((3,1),dtype='int')

print('Matrix L:')
print(L)
print('Matrix M:')
print(M)
print('Matrix N:')
print(N)

Output:

Matrix L:
[0. 0. 0.]
Matrix M:
[[0. 0. 0. 0.]
 [0. 0. 0. 0.]]
Matrix N:
[[0]
 [0]
 [0]]

Here the null matrix $N$ is also a Column Matrix.

There is one more function to create a zero matrix, with the same number of rows and columns as any given matrix called numpy.zero_like( ). The same data type as the input matrix is followed.

A = np.array([[1,2,3],[4,5,6]])
K = np.zeros_like(A)

print('Matrix A:')
print(A)
print('Zero Matrix K:')
print(K)

Output:

Matrix A:
[[1 2 3]
 [4 5 6]]
Zero Matrix K:
[[0 0 0]
 [0 0 0]]

Matrix of Ones

The numpy.ones( ) function is used to create a matrix made up of ones. Similar to the zeros function, this function takes shape of a matrix as an input and the data type is an optional argument.

X = np.ones(5)
Y = np.ones((3,2), dtype='complex')
Z = np.ones((4,1), dtype='bool')

print('Matrix X:')
print(X)
print('Matrix Y:')
print(Y)
print('Matrix Z:')
print(Z)

Output:

Matrix X:
[1. 1. 1. 1. 1.]
Matrix Y:
[[1.+0.j 1.+0.j]
 [1.+0.j 1.+0.j]
 [1.+0.j 1.+0.j]]
Matrix Z:
[[ True]
 [ True]
 [ True]
 [ True]]

As you can see from the $Z$ matrix, we can create a matrix with the boolean value True as elements. Also, we can create a matrix with boolean value False using the numpy.zeros( ) function, as 0 corresponds to boolean False and 1 corresponds to boolean True.

This function also has the numpy.ones_like( ) variant to create a matrix of ones with the same order and data type as the given matrix.

Identity Matrix

There are two functions available in the NumPy library to create an Identity Matrix, which has diagonal elements as ones and the remaining elements as zeros.

The first way is to use the numpy.identity( ) function. It takes an integer as an argument, which specifies the order of the identity matrix. The data type (dtype) argument is optional.

I3 = np.identity(3)
I4 = np.identity(4,dtype='int')

print('Identity Matrix of Order 3:')
print(I3)
print('Identity Matrix of Order 4:')
print(I4)

Output:

Identity Matrix of Order 3:
[[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]
Identity Matrix of Order 4:
[[1 0 0 0]
 [0 1 0 0]
 [0 0 1 0]
 [0 0 0 1]]

We can only create the standard square identity matrix with the numpy.identity( ) function. A more versatile, numpy.eye( ) function is also available with few more arguments.

The numpy.eye( ) function takes main four arguments, three of which are optional. The first argument is the number of rows ($N$).

If the second argument ($M$), i.e., the number of columns is specified, then we get a $N \times M$ rectangular identity-like matrix. Otherwise, we get a standard identity matrix of the order $N$.

The numpy.eye( ) function takes two more optional arguments. The first is $k$, which specifies the offset of the diagonal. The default value of $k$ is 0. The positive value of $k$ shifts the non-zero diagonal above the principal diagonal and the negative value shifts it below.

The next optional argument is the data type with float as a default. Let’s see a few examples:

I3 = np.eye(3)
P = np.eye(3,4,dtype='int')
Q = np.eye(5,k=1)
R = np.eye(4,k=-1)

print('Identity Matrix of Order 3:')
print(I3)
print('Matrix P:')
print(P)
print('Matrix Q:')
print(Q)
print('Matrix R:')
print(R)

Output:

Identity Matrix of Order 3:
[[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]
Matrix P:
[[1 0 0 0]
 [0 1 0 0]
 [0 0 1 0]]
Matrix Q:
[[0. 1. 0. 0. 0.]
 [0. 0. 1. 0. 0.]
 [0. 0. 0. 1. 0.]
 [0. 0. 0. 0. 1.]
 [0. 0. 0. 0. 0.]]
Matrix R:
[[0. 0. 0. 0.]
 [1. 0. 0. 0.]
 [0. 1. 0. 0.]
 [0. 0. 1. 0.]]

The identity matrix is really special because it satisfies the conditions for multiple types of matrices.

Creating different types of matrices

Matrix full of a given number

Similar to the matrix of ones, we can create a matrix with all the elements having the same desired value. For this purpose, we can make use of the numpy.full( ) function. This function requires the shape and the fill_value to be specified. The dtype is optional.

X = np.full((2,3),8)
Y = np.full((2,2),np.inf)

print('Matrix X:')
print(X) 
print('The infinity matrix of order 2:')
print(Y)

Output:

Matrix X:
[[8 8 8]
 [8 8 8]]
The infinity matrix of order 2:
[[inf inf]
 [inf inf]]

Similar to the numpy.zeros_like( ) and numpy.ones_like( ), we have the numpy.full_like( ) function available to create a matrix full of given number and having the same shape and data type as any given matrix.

Diagonal Matrix

The Diagonal Matrix has the non-zero elements along the principal diagonal and all the other elements are zeros. We can create this type of matrix using the numpy.diag( ) function, which takes the diagonal elements as the input. This function also takes an optional argument $k$, which is used to select the diagonal above ($k$ > 0) or below ($k$ < 0) the principal diagonal to be filled. The default value of $k$ is zero indicating the principal diagonal.

Just to remind you that the diagonal starting from the top-left corner is the principal diagonal.

my_list = [4,9,6]

L = np.diag(my_list)
M = np.diag(my_list,k=1)
N = np.diag(my_list,k=-1)

print('Matrix L:')
print(L)
print('Matrix M:')
print(M)
print('Matrix N:') 
print(N)

Output:

Matrix L:
[[4 0 0]
 [0 9 0]
 [0 0 6]]
Matrix M:
[[0 4 0 0]
 [0 0 9 0]
 [0 0 0 6]
 [0 0 0 0]]
Matrix N:
[[0 0 0 0]
 [4 0 0 0]
 [0 9 0 0]
 [0 0 6 0]]

And interestingly enough, we can use the same function to extract the diagonal elements of any given matrix. We can specify which diagonal to extract with the optional argument $k$. By default, it extracts the principal diagonal.

Note that the output is also a NumPy array (NumPy ndarray object). Let’s see a few examples of that.

A = np.array([[2,3,5],[7,11,13],[17,19,23]])
B = np.diag(A)
C = np.diag(A,k=1)
D = np.diag(A,k=-1)

print('Matrix A:')
print(A)
print('Principal diagonal of A:')
print(B)
print('Diagonal above the principal diagonal of A:')
print(C)
print('Diagonal below the principal diagonal of A:')
print(D)

Output:

Matrix A:
[[ 2  3  5]
 [ 7 11 13]
 [17 19 23]]
Principal diagonal of A:
[ 2 11 23]
Diagonal above the principal diagonal of A:
[ 3 13]
Diagonal below the principal diagonal of A:
[ 7 19]

There is one more function called numpy.diagflat(), which serves the same purpose as the numpy.diag( ) function. The only difference is that we can use a 2D array (Rectangular or Square Matrix) as the input and the function flattens it to be a 1D array [Row Matrix] to get the diagonal elements.

X = np.array([[9,8],[7,6]])

Y = np.diagflat(X)
Z = np.diagflat(X,1)

print('Input Matrix X:')
print(X)
print('Input Matrix Y:')
print(Y)
print('Input Matrix Z:')
print(Z)

Output:

Input Matrix X:
[[9 8]
 [7 6]]
Input Matrix Y:
[[9 0 0 0]
 [0 8 0 0]
 [0 0 7 0]
 [0 0 0 6]]
Input Matrix Z:
[[0 9 0 0 0]
 [0 0 8 0 0]
 [0 0 0 7 0]
 [0 0 0 0 6]
 [0 0 0 0 0]]

Triangular Matrix

As per the standard definition, the square matrix is said to be a Triangular Matrix if the elements on one side of the principal diagonal are zeros.

The NumPy library offers three functions to create the triangular matrices namely numpy.tri( ), numpy.tril( ) and numpy.triu( ). The numpy.tri( ) function creates a matrix with the elements at and below the selected diagonal as ones and the remaining elements as zeros.

It takes the number of rows ($N$) and the number of columns ($M$) as the input arguments ($M$ = $N$ if not specified). The diagonal is selected using the optional argument $k$.

Similar to the other functions, $k$ = 0 is the default setting and it selects the principal diagonal while $k$ > 0 and $k$ < 0 select the diagonals above and below the principal diagonal, respectively. The function also has an optional dtype argument.

P = np.tri(3)
Q = np.tri(4,3,dtype=int)
R = np.tri(3,k=1)

print('Matrix P:')
print(P)
print('Matrix Q:')
print(Q)
print('Matrix R:')
print(R)

Output:

Matrix P:
[[1. 0. 0.]
 [1. 1. 0.]
 [1. 1. 1.]]
Matrix Q:
[[1 0 0]
 [1 1 0]
 [1 1 1]
 [1 1 1]]
Matrix R:
[[1. 1. 0.]
 [1. 1. 1.]
 [1. 1. 1.]]

The numpy.tril( ) function returns the Lower Triangular Matrix while the numpy.triu( ) function returns the Upper Triangular Matrix of the given input matrix.

The optional argument $k$ works the same way as before. When a diagonal is selected using $k$, the numpy.tril( ) function returns the copy of the input matrix with the elements above the selected diagonal zeroed.

On the other hand, the numpy.triu( ) function returns the copy of the input matrix with the elements below the $k^{th}$ diagonal zeroed.

A = np.reshape(np.array(range(1,10)),(3,3))
L = np.tril(A)
U = np.triu(A)

P = np.tril(A,k=-1)
Q = np.triu(A,k=-1)

print('Matrix A:')
print(A)
print('------------------------------')
print('Lower Triangular Matrix L:')
print(L)
print('Upper Triangular Matrix U:')
print(U)
print('------------------------------')
print('Matrix P:')
print(P)
print('Matrix Q:')
print(Q)

Output:

Matrix A:
[[1 2 3]
 [4 5 6]
 [7 8 9]]
------------------------------
Lower Triangular Matrix L:
[[1 0 0]
 [4 5 0]
 [7 8 9]]
Upper Triangular Matrix U:
[[1 2 3]
 [0 5 6]
 [0 0 9]]
------------------------------
Matrix P:
[[0 0 0]
 [4 0 0]
 [7 8 0]]
Matrix Q:
[[1 2 3]
 [4 5 6]
 [0 8 9]]

Random Matrix

We can create matrices with the random numbers using the numpy.random module. The numpy.random.rand( ) function from this module generates the matrix with the pseudo-random elements with the uniform distribution over 0 (inclusive) to 1 (exclusive).

This function takes a tuple of the desired output dimension as an input.

R1 = np.random.rand(3)
R2 = np.random.rand(3,2)

print('Random Row Matrix R1:')
print(R1)
print('------------------------------')
print('Random Matrix R2:')
print(R2)

Output:

Random Row Matrix R1:
[0.54874482 0.94463367 0.61183853]
------------------------------
Random Matrix R2:
[[0.19483878 0.99341358]
 [0.7436171  0.2471102 ]
 [0.13959398 0.56015437]]

The numpy.random.randint( ) function generates the matrix with the pseudo-random integer elements within the low (inclusive) to high (exclusive) bounds. If high is not specified, the function uses the range [0, low).

RI1 = np.random.randint(1,7,size=3)
RI2 = np.random.randint(10,size=(3,2))
RI3 = np.random.randint(10,100,size=(4,3))

print('Three dice are rolled (RI1):' )
print(RI1)
print('Random Matrix RI2:')
print(RI2)    
print('Random Matrix RI3:')
print(RI3)

Output:

Three dice are rolled (RI1):
[4 1 1]
Random Matrix RI2:
[[2 4]
 [9 9]
 [4 7]]
Random Matrix RI3:
[[97 14 45]
 [28 56 26]
 [16 41 21]
 [20 20 96]]

You can check the articles in my Linear Algebra Using Python series to see the uses of the above matrix creation routines. You can also check the Github repository

Please let me know if you find the article useful. Please comment with questions or suggestions and share the article with your friends.