Monthly Archives for October 2020

Solving symbolic matrix equations in Python with SymPy

I am no Python expert, and have only recently encountered SymPy, for symbolic calculations. I have been trying to do some (relatively simple) matrix calculations, and it has taking me an embarrassingly long time to figure out how to do this. So, I am sharing what I have learned here to help someone else avoid the rather large number of internet searches I had to do to piece it together.

There is a lot out there on how to use SymPy to solve matrix equations of the form Ax = b. What I am interested in is taking a bunch of given matrices (with numerical values) and constants, performing some operations with an unknown matrix, and setting each entry of the final matrix equal to zero and solving. In other words, suppose we are given n \times n matrices A and B, which are determined beforehand. There is some constant t which can be varied (this forms part of an iterative scheme), and an unknown matrix C, which is represented purely symbolically, as such:

\begin{bmatrix} c_{0,0} & c_{0,1} & \cdots &c_{0, n-1} \\ \vdots & \vdots & \ddots & \vdots \\ c_{n-1,0} & c_{n-1,1} & \cdots & c_{n-1,n-1} \end{bmatrix}

There is a function f_{A,B,t}: \mathbb{R}^{n\times n}\to\mathbb{R}^{n\times n}, and we want to find the entries of C for which

f_{A,B,t}(C) = \mathbf{0},

where \mathbf{0} indicates the n\times n zero matrix.

In order to solve an equation in SymPy, you have to declare the “symbols” that you are solving for. Now, defining a matrix symbol in SymPy is easy, but this did not help me in solving for the matrix, and I kept getting an empty output. I needed a way to iteratively declare each entry of the matrix as a symbol, whilst putting them together as a single matrix. This turned out to be the key to the whole thing. First, let us state the preamble:

import sympy as sym
from sympy.interactive import printing
printing.init_printing(use_latex=True)

The last two lines have no effect on the calculations, but they do give the option of displaying your matrices very nicely. Next, let us define some function with which to work:

def matrixfunction(A,B,C,t):
    mf = A - t*C + B.inv()
    return mf

(The final part of the last line is simply how we compute the inverse of B.) Here is one of the first things that tripped me up. In SymPy, you should distinguish between operations involving symbolic matrices and usual operations between matrices. For instance, if I were to declare D and E to be two arbitrary 5\times 5 matrices and wanted, for instance, to multiply them, I would use

D = sym.MatrixSymbol('D', 5, 5)
E = sym.MatrixSymbol('E', 5, 5)
sym.MatMul(D,E)

The output to this would be 

D*E

and we would be able to see the symbolic entries of this matrix by using

X = sym.MatMul(D,E)
X.as_explicit()

The same holds for MatAdd. However, if you have defined the matrix by declaring all of its entries to be symbols, there does not seem to be a need to use this method, and a simple * can be used for multiplication and + for addition.

Our objective is now to set each entry in the matrix obtained from the function “matrixfunction” equal to zero and solve for the unknown matrix C. To do so, we define 

def solvefor(A,B,t):
    C = Matrix(n,n,sym.symbols('D0:n(0:n)'))
    sol = sym.solve(matrixfunction(A,B,C,t))
    display sol

Of course, in the above, n needs to be replaced by an actual numerical value. Once the functions are defined, we can assign values to A, B and t and run “solvefor”. The output will be a set of values assigning the solution value to each entry of C. In this case, B has to be invertible. Note that SymPy automatically sets the argument of “sym.solve” equal to zero unless otherwise instructed – in this case, it is set equal to the zero matrix. We are also not specifying which symbols to solve for, since SymPy will automatically solve for the free variables here.