An Introduction to Matrix Computations in Python: Just the basics

We can compute with matrices in both the Sympy and Numpy libraries. Numerical Python is intended for numerical computations and estimates. Symbolic Python can do symbolic computations like rref() or compute with matrices having undetermined variables in the mix. NumPy automatically treats matrices like a 2D array, that is, a higher dimensional vector. SymPy defaults to working with matrices from the linear transformation perspective.

After some internal debate, I decided to intermingle NumPy and SymPy syntax for working with matrices. Please give feedback if it seems this section should be re-organized into NumPy and Matrics followed by SymPy and Matrices.

Some matrices are frequently needed for computations. Common matrices that we can generate automatically include a (multiplicative) identity matrix, a zero matrix, and a matrix of all ones.

The commands for the \(n\times n\) multiplicative identity matrix are different: SymPy’s eye(n) vs NumPy’s identity(n).

The commands for matrices of ones() and zeros() are identical, and take in the number of rows and columns as arguments. For NumPy, however, you must enter the shape of the array as an ordered pair with parentheses.

You can add, subtract and scale matrices using both NumPy and SymPy.

Recall that matrix multiplication is more sophisticated than the Hadamard product, a basic element-wise multiplication of two arrays, since it developed to be independent of coordinates and correspond to function composition. Recall further that NumPy automatically treats matrices like a 2D array, that is, a higher dimensional vector. If you naively multiply NumPy arrays you get the Hadamard product, rather than traditional matrix multiplication (try below).

If instead you want matrix multiplication, you tell NumPy to multiply using the matmul() command, or as a short-hand using @ (try below).

Meanwhile SymPy defaults to working with matrices from the linear transformation perspective giving you matrix multiplication by default and requiring us to specify any Hadamard product, or element-wise multiplication.

Square matrices can be multiplied by themselves and at times we want to multiply repeatedly, e.g. \(A^3=A\cdot A \cdot A\text{.}\) In SymPy we can perform repeated matrix multiplication using Python’s exponent operator A**3. In NumPy we can call linalg.matrixpower(matrix, integer power).

Recall that we can find multiplicative matrix inverses by hand by augmenting a matrix with the identity matrix and row reducing to the identity. Thus we can find a matrix inverse in SymPy simply by applying rref() to a matrix augmented with the identity matrix, which we cannot do with NumPy.

Both SymPy and NumPy have a dedicated matrix inverse command, however, NumPy’s linalg.inv() and SymPy’s inv().

We can solve matrix equations using matrix inverses.

While inverting a matrix helps us solve matrix equations symbolically and develop theory, in practice, we seldom invert a matrix with technology in order to solve systems of linear equations. Some matrices are ill-conditioned in that small changes what we multiply them by can create big changes in the product. Since rounding produces small changes, ill-conditioned matrices can create terminal errors in our solutions.

In general, matrices that are nearly singular, that is, almost not invertible, can be ill-conditioned. In Section 2.4 we will learn how to use determinants to discover whether or not a matrix is singular or close to singular.