Lesson 6. Matrices

2023-01-21

Numerical Linear Algebra

Computers enable the solution of large scale problems in linear algebra. A modern desktop computer can solve a \(5,000\times5,000\) linear system in under half a second, and find all of the eigenvalues and eigenvectors of a \(5,000\times5,000\) real symmetric matrix in about 12 seconds. These and similar capabilities have opened up many new application areas for mathematical modeling and analysis.

Objectives

In this lesson you will learn how to work with matrices in Python. When you complete the lesson, you should know

how to create matrices;
convenient ways to construct commonly-occuring kinds of matrices;
indexing techniques for elements, rows, columns and sub-matrices;
how to perform matrix operations in Python;
how to solve square linear systems.

Creating Matrices

In Python, the best way to represent a matrix is as a two-dimensional NumPy array. For small matrices, you can create the associated array by forming a list of the rows of the matrix, with each row consisting of a list of numbers. For example, the \(3\times4\) matrix \[ A=\begin{array}{cccc} 2& 0&-7& 1\\ 5& 9& 3&-1\\ 4& 1& 0&-6\end{array} \] can be constructed by forming the list (of lists)

l = [ [2, 0, -7,  1], 
      [5, 9,  3, -1],
      [4, 1,  0, -6] ]

and doing

A = array(l)

The following NumPy functions are often useful when solving matrix problems:

zeros(m, n) returns an \(m\times n\) array of zeros;
ones(m, n) returns and \(m\times n\) array of ones;
eye(n) return an \(n\times n\) matrix with ones on the main diagonal, and zeros everywhere else (i.e., an identity matrix).

The module numpy.random provides functions for constructing random matrices:

rand(m, n) returns an \(m\times n\) matrix whose entries are uniformly distributed in the unit interval \([0,1]\);
randn(m, n) returns an \(m\times n\) matrix whose entries are normally distributed with mean 0 and variance 1.

Indexing

To reference an entry in a two dimensional array A we use two indices enclosed in square brackets, generalising the syntax used for a one dimensional array. The language of matrices is commonly employed, so that the entry A[i,j] is said to have row index i and column index j, and we visualise the entries arranged as in a matrix. Although we commonly use the terms ‘matrix’ and ‘two-dimensional array’ interchangeably, note that, unlike the usual practice in mathematics, Python indices start at 0 instead of 1. Thus, the entries of an \(m\times n\) array are arranged as follows \[ \begin{array}{cccc} A[0,0]&A[0,1]&\cdots&A[0,n-1]\\ A[1,0]&A[1,1]&\cdots&A[1,n-1]\\ \vdots&\vdots&\ddots&\vdots\\ A[m-1,0]&A[m-1,1]&\cdots&A[m-1,n-1]. \end{array} \]

Sub-arrays of a matrix can be defined by slicing. The ith row of A is referenced by A[i,:], and the jth column by A[:,j]. In general, A[i1:i2,j1:j2] is the sub-matrix with entries A[i,j] for i from i1 up to, but not including, i2, and for j from j1 up to, but not including, j2. Thus, for the \(3\times4\) array considered earlier, \[ A=\begin{array}{cccc} 2& 0&-7& 1\\ 5& 9& 3&-1\\ 4& 1& 0&-6\end{array} \] we have \[ A[1,:]=\begin{array}{cccc} 5&9&3&-1 \end{array},\qquad A[:,2]=\begin{array}{c} -7\\ 3\\ 0 \end{array},\qquad A[1:3,2:4]=\begin{array}{cc} 3&-1\\ 0&-6\end{array}. \]

The ith row of a matrix A can be referenced as A[i] as an alternative to A[i,:]. This syntax reflects the idea that A is defined from a list of rows.

A sub-array is itself an array with entries indexed from \(0\). For example, if A is as above and we put

B = A[1:3,2:4]

then B[0,0]=3, B[0,1]=-1, B[1,0]=0 and B[1,1]=-6. Moreover, in this situation B references the same memory locations as the corresponding entries of A, so that, e.g., if we do

B[1,1] = 8

then the value of A[2,3] is also changed from -6 to 8.

The dimensions of any NumPy array A can be extracted from the tuple A.shape. Thus, with A defined as above, A.shape = (3,4). For a general 2-dimensional array A, the number of rows equals A.shape[0] and the number of columns equals A.shape[1].

Matrix Arithmetic

Addition of two-dimensional arrays is defined elementwise, consistent with the way we define addition of matrices in linear algebra. Likewise for the product of a scalar and a two-dimensional array. Thus, if

A = array([[1, 2],
           [3, 4],
           [5, 6]])
B = -2 * A
C = 2 * A + B

then C will be the \(3\times2\) zero matrix.

However, matrix multiplication differs from elementwise multiplication of two-dimensional arrays. NumPy uses @ for the former operation and * the latter. For example, if

A = array([[3, -2],
           [5,  1]])
B = array([[1,  0],
           [7, -4]])
C = A @ B
D = A * B

then \[ C = \begin{array}{cc} -11 & 8\\ 12 & -4\end{array},\qquad D = \begin{array}{cc} 3& 0\\ 35& -4\end{array}. \]

In general, A @ B is defined iff the number of columns in the matrix A equals the number of rows in the matrix B, or in other words, iff A.shape[1] equals B.shape[0]. The pointwise product A * B is defined iff both arrays have the same dimensions, or in other words iff A.shape equals B.shape.

Matrices and Vectors

Type the following commands in the IPython pane.

A = array([[3, -2],
           [1,  5]])
x = array([1, 2])
x_col = array([[1], 
               [2]])
y = A @ x
y_col = A @ x_col

Notice that y is a one-dimensional array whereas y_col is a \(2\times1\) Matrix, and thus a column vector.

Exercise. Compare y[0] and y_col[0].

Exercise. Define the row vector x_row = array([[1, 2]]) and verify that A @ x_row is undefined but x_row @ A is the mathematically correct row vector.

For any matrix A, NumPy uses the syntax A.T for the transpose of A. Thus, if A is as in the preceding example, then A.T is \[ \begin{array}{cc} 3 & 1 \\ -2 & 5\end{array}. \] Similarly, the column vector x_col equals x_row.T, and the row vector x_row equals x_col.T.

Linear Systems

Consider the following \(3\times3\) system of linear equations, \[ \begin{array}{cc} 5x_1+ x_2+2x_3&=17,\\ -2x_1+9x_2+4x_3&=-11,\\ 3x_1+ x_2+8x_3&=43. \end{array} \] To compute the solution, we define the coefficient matrix A and right-hand side vector b,

A = array([[ 5,  1,  2],
           [-2,  9,  4],
           [ 3,  1,  8]])
b = array([[ 17],
           [-11],
           [ 43]])

and then use the solve function from the linalg submodule of SciPy.

from scipy.linalg import solve
x = solve(A, b)
r = b - A @ x

Here, we compute the residual vector \(r=b-Ax\) as a simple check that the solve function worked as expected. In this case, you will find that r is the zero vector, but in general you have to expect the entries of r to be non-zero but very small.

Calling solve with a singular matrix raises LinAlgError: Matrix is singular.

Summary

This lesson taught you how to

create an array from a list of lists;
use zeros or ones to create a matrix of all zeros or all ones;
reference a row or column of a matrix by slicing;
assign a variable to a sub-matrix;
perform matrix arithmetic;
use scipy.linalg.solve to solve a square, non-singular system of linear equations.