Solving Linear Systems (SymPy)

Most linear algebra courses start by considering how to solve a system of linear equations.

$$\begin{array}{rl}{a}_{0,0}{x}_0+a_{0,1}{x}_0+\cdots a_{0,n-1}{x}_0& =b_0\\ a_{1,0}{x}_0+a_{1,1}{x}_0+\cdots a_{1,n-1}{x}_0& =b_1\\ & ⋮\\ a_{n-1,0}{x}_0+a_{n-1,1}{x}_0+\cdots a_{n-1,n-1}{x}_0& =b_{m-1}\end{array}$$

In matrix notation, this linear system of $n$ equations in $n$ unknowns ($x_0,\dots ,x_{n-1}$), can be written as

$$A\overrightarrow{x}=\overrightarrow{b}$$

where

$$A=\left[\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}\\ ⋮& ⋮& & ⋮\\ a_{n-1,0}& a_{n-1,1}& \cdots & a_{n-1,n-1}\end{array}\right],\overrightarrow{x}=\left[\begin{array}{c}{x}_0\\ x_1\\ ⋮\\ x_{n-1}\end{array}\right],\overrightarrow{b}=\left[\begin{array}{c}{b}_0\\ b_1\\ ⋮\\ b_{n-1}\end{array}\right]$$

We'll look at how to use SymPy to solve linear systems, including those with free variables. We'll also look at how to do row operations (Gaussian elimination).

First we import the necessary libraries.

from sympy import *
from sympy.solvers.solveset import linsolve
x, y, z = symbols('x, y, z')

.solve and .linsolve

Let's start with a 3x3 linear system:

$$\begin{array}{rl}{x}_1+x_2-2x_3& =4\\ 2x_2+x_3& =3\\ 2x_1+x_2-5x_3& =7\end{array}$$

Verify (by substitution) that $\overrightarrow{x}=[0,2,-1]$ is a solution.

System of Equations

To solve a linear system represented by a system of linear equations we use the sympy.solve funtion.

eq1 = x + y - 2*z
eq2 = 2*y + z
eq3 = 2*x + y - 5*z
solve([eq1-4, eq2-3, eq3-7], (x, y, z)) # (1)

1. $A$ may can be a rectangular matrix.

{x: 0, y: 2, z: -1}

Matrix Form

For a linear system expressed in matrix form $A\overrightarrow{x}=\overrightarrow{b}$ we can use linsolve.

M = Matrix([[1,1,-2],[0,2,1],[2,1,-5]])
b = Matrix([4,3,7])
linsolve([M,b], x, y, z)

{ (0,2,-1) }

Augmented Matrix Form

For a linear system expressed as an augmented matrix $[A|\overrightarrow{b}]$ we can use linsolve.

M = Matrix([[1,1,-2],[0,2,1],[2,1,-5]])
b = Matrix([4,3,7])
M_aug = M.row_join(b)
linsolve(M_aug, (x, y, z))

{ (0,2,-1) }

System with Free Variables

The three methods above can be used to solve a system with free variables.

Consider the system

$$\begin{array}{rl}x+y+z& =1\\ x+y+2z& =3\end{array}$$

which has solution space $\{(-t-1,t,2):t\in \mathbb{R}\}$.

equations$Ax=b$augmented matrix

eq1 = x + y + z
eq2 = x + y + 2*z
solve([eq1-1, eq2-3], (x, y, z))

{x: -y - 1, z: 2}

M = Matrix([[1,1,1],[1,1,2]])
b = Matrix([1,3])
linsolve([M,b], (x, y, z))

{ ( -y-1, y, 2) }

M = Matrix([[1,1,1],[1,1,2]])
b = Matrix([1,3])
M_aug = M.row_join(b)
linsolve(M_aug, (x, y, z))

{ ( -y-1, y, 2) }

Example

Consider another example: $$2x_1+3x_2=11$$ $$3x_1+6x_2=7$$

A = Matrix([[2,3],[3,6]])
b = Matrix([11,7])
linsolve(A,b,(x,y))

$\left\{(15,-\frac{19}{3})\right\}$

We can verify this answer is correct by observing $A\overrightarrow{x}=\overrightarrow{b}$.

sol=linsolve([A,b],(x,y))
for s in sol:
    pprint(A*Matrix(s)) # (1)

1. pprint is pretty print, which displays matrices in the proper format.

$\left[\begin{array}{c}11\\ 7\end{array}\right]$

Gaussian Elimination

The general procedure learned to solve a system of linear equations is Gaussian elimination. The goal is to apply row operations to "eliminate" all the variables except for one in each equation (if possible), and this resulting linear system has the same set of solutions as the original, but is easier to solve. We generally do this on the matrix form of the system, and call this the reduced row echelon form (rref) of the matrix. The elementary row operations we use are:

Elementary Row Operations

1. Replace row $j$ with $k$ times row $i$ added to row $j$. $R_j\leftarrow k{R}_i+R_j$
2. Multiply row $i$ by a scalar $k$. $R_i\leftarrow k{R}_i$
3. Swap rows $i$ and $j$. $R_i\leftrightarrow R_j$

For example, consider the matrix

$$A=\left[\begin{array}{ccc}2& 1& 0\\ 1& 1& -1\\ 0& -1& 2\end{array}\right]$$

Since we are working in python we will index rows starting at $0$. In other words, an $n\times m$ matrix has rows $R_0$ to $R_{n-1}$.

Here are the row operations taking it to rref:

$$A=\left[\begin{array}{ccc}2& 1& 0\\ 1& 1& -1\\ 0& -1& 2\end{array}\right]\stackrel{R_0\leftrightarrow R_1}{\to }\left[\begin{array}{ccc}1& 1& -1\\ 2& 1& 0\\ 0& -1& 2\end{array}\right]\stackrel{R_1\leftarrow -2R_0+R_1}{\to }\left[\begin{array}{ccc}1& 1& -1\\ 0& -1& 2\\ 0& -1& 2\end{array}\right]$$

$$\stackrel{R_2\leftarrow -R_1+R_2}{\to }\left[\begin{array}{ccc}1& 1& -1\\ 0& -1& 2\\ 0& 0& 0\end{array}\right]\stackrel{R_0\leftarrow R_1+R_0}{\to }\left[\begin{array}{ccc}1& 0& 1\\ 0& -1& 2\\ 0& 0& 0\end{array}\right]\stackrel{R_1\leftarrow -R_1}{\to }\left[\begin{array}{ccc}1& 0& 1\\ 0& 1& -2\\ 0& 0& 0\end{array}\right]$$

To do row operations in SymPy us the elemetary_row_op method on a Matrix object.

elementary_row_op(op='n->kn', row=None, k=None, row1=None, row2=None)

Performs the elementary row operation

• "n->kn" (row n goes to k*n)
• "n<->m" (swap row n and row m)
• "n->n+km" (row n goes to row n + k*row m)
• row is the index of the row to be replaced.
• row1 is the row that will be added to row
• row1 and row2 are the two rows in the swap.

Example

We'll redo the example above, this time using SymPy.

A = Matrix([[2,1,0],[1,1,-1],[0,-1,2]])
pprint(A)
print("Swap rows 0 and 1:")
A1 = A.elementary_row_op("n<->m",0,1)
pprint(A1)
print("Add -2 times row 0 to row 1:")
A2 = A1.elementary_row_op("n->n+km",row=1,k=-2,row1=0)
pprint(A2)
print("Add -1 times row 1 to row 2:")
A3 = A2.elementary_row_op("n->n+km",row=2,k=-1,row1=1)
pprint(A3)
print("Add 1 times row 1 to row 0:")
A4 = A3.elementary_row_op("n->n+km",row=0,k=1,row1=1)
pprint(A4)
print("Multiply row 1 by -1:")
A5 = A4.elementary_row_op("n->kn",row=1,k=-1)
pprint(A5)

 [ 2  1  0]
 [ 1  1 -1]
 [ 0 -1  2]
Swap rows 0 and 1:
 [ 1  1 -1 ]
 [ 2  1  0 ]
 [ 0 -1  2 ]
Add -2 times row 0 to row 1:
 [ 1  1 -1 ]
 [ 0 -1  2 ]
 [ 0 -1  2 ]
Add -1 times row 1 to row 2:
 [ 1  1 -1 ]
 [ 0 -1  2 ]
 [ 0  0  0 ]
Add 1 times row 1 to row 0:
 [ 1  0  1 ]
 [ 0 -1  2 ]
 [ 0  0  0 ]
Multiply row 1 by -1:
 [ 1  0  1 ]
 [ 0  1 -2 ]
 [ 0  0  0 ]

Solving a linear system with Gaussian Elimination

We'll consider the following linear system which has more unknowns than equations. This means the coefficient matrix will not be square.

$$\begin{array}{rl}2x_1+4x_3+4x_4& =14\\ -x_1-x_3& =-3\\ -3x_1+6x_4& =3\end{array}$$

Write it as an augmented matrix:

$$\left[\begin{array}{ccccc}2& 0& 4& 4& 14\\ -1& 0& -1& 0& -3\\ -3& 0& 0& 6& 3\end{array}\right]$$

Use row operations to take it to rref.

A = Matrix([[2,0,4,4],[-1,0,-1,0],[-3,0,0,6]])
b = Matrix([14,-3,3])
aug_mat = A.row_join(b) # augmented matrix
aug_mat

$\left[\begin{array}{ccccc}2& 0& 4& 4& 14\\ -1& 0& -1& 0& -3\\ -3& 0& 0& 6& 3\end{array}\right]$

A1 = aug_mat.elementary_row_op("n<->m",0,1)
A1

$\left[\begin{array}{ccccc}-1& 0& -1& 0& -3\\ 2& 0& 4& 4& 14\\ -3& 0& 0& 6& 3\end{array}\right]$

A2 = A1.elementary_row_op("n->n+km",row=1,k=2,row1=0)
A3 = A2.elementary_row_op("n->n+km",row=2,k=-3,row1=0)
A4 = A3.elementary_row_op("n->kn",row=0,k=-1)
A4

$\left[\begin{array}{ccccc}1& 0& 1& 0& 3\\ 0& 0& 2& 4& 8\\ 0& 0& 3& 6& 12\end{array}\right]$

A5 = A4.elementary_row_op("n->kn",row=2,k=1/3) # (1)
A5

1. Use k=Rational(1,3) to avoid using floating point numbers.

$\left[\begin{array}{ccccc}1& 0& 1& 0& 3\\ 0& 0& 2& 4& 8\\ 0& 0& 1& 2& 4\end{array}\right]$

A6 = A5.elementary_row_op("n->n+km",row=1,k=-2,row1=2)
A7 = A6.elementary_row_op("n->n+km",row=0,k=-1,row1=2)
A7

$\left[\begin{array}{ccccc}1& 0& 0& -2& -1\\ 0& 0& 0& 0& 0\\ 0& 0& 1& 2& 4\end{array}\right]$

A8 = A7.elementary_row_op("n<->m", 1,2)
A8

$\left[\begin{array}{ccccc}1& 0& 0& -2& -1\\ 0& 0& 1& 2& 4\\ 0& 0& 0& 0& 0\end{array}\right]$

At this point we have reduced the matrix to rref. We see there are two free variables ($x_2$ and $x_4$). The general solution is:

$$\overrightarrow{x}=\left[\begin{array}{c}-1\\ 0\\ 4\\ 0\end{array}\right]+s\left[\begin{array}{c}0\\ 1\\ 0\\ 0\end{array}\right]+t\left[\begin{array}{c}2\\ 0\\ -2\\ 1\end{array}\right],\text{ for }s,t\in \mathbb{R}$$

.rref

SymPy has a built-in function to put a matrix in reduced row echelon form (rref). To put a matrix into reduced row echelon form use the method .rref. .rref returns a tuple of two elements. The first is the reduced row echelon form, and the second is a tuple of indices of the pivot columns. The length of this tuple is the rank of the matrix.

Example

Compute the rref of the matrix

$$\left[\begin{array}{ccc}2& 1& 0\\ 1& 1& -1\\ 0& -1& 2\end{array}\right]$$

A = Matrix([[2,1,0],[1,1,-1],[0,-1,2]])
pprint(A.rref())

$(\left[\begin{array}{ccc}1& 0& 1\\ 0& 1& -2\\ 0& 0& 0\end{array}\right],(0,1))$

Example

Compute the rref of the matrix

$$\left[\begin{array}{cccc}2& 0& 4& 4\\ -1& 0& -1& 0\\ -3& 0& 0& 6\end{array}\right]$$

A = Matrix([[2,0,4,4],[-1,0,-1,0],[-3,0,0,6]])
pprint(A.rref())

$(\left[\begin{array}{cccc}1& 0& 0& -2\\ 0& 0& 1& 2\\ 0& 0& 0& 0\end{array}\right],(0,2))$

Example

For the linear system, find the rref of the augmented matrix.

$$\left[\begin{array}{cccc}2& 0& 4& 4\\ -1& 0& -1& 0\\ -3& 0& 0& 6\end{array}\right]\overrightarrow{x}=\left[\begin{array}{c}14\\ -3\\ 3\end{array}\right]$$

A = Matrix([[2,0,4,4],[-1,0,-1,0],[-3,0,0,6]])
b = Matrix([14,-3,3])
aug_mat = A.row_join(b) # augmented matrix
pprint(aug_mat.rref())

$(\left[\begin{array}{ccccc}1& 0& 0& -2& -1\\ 0& 0& 1& 2& 4\\ 0& 0& 0& 0& 0\end{array}\right],(0,2))$

Computing an Inverse with Gaussian Elimination

Recall, if we have an invertible matrix $A$, then we can compute the inverse using Gaussian elimination as follows:

$$[A|I]\underset{\text{to rref}}{\overset{\text{gaussian elim}}{\to }}[I|B]$$

Matrix $B$ is the inverse of $A$.

Note

If you are wondering why the inverse matrix can be computed in this way, it comes back to elementary matrices. Suppose $n$ row operations are applied to take $A$ to the identity $I$. Let $E_1,E_2,\dots E_n$ be the $n$ elementary matrices corresponding to these row operations. Then, since row operations could be done by multiplying by elementary matrices, it follows that

$$E_n\cdots E_2{E}_{1}A=I$$

Therefore, multiplying $A^{-1}$ to the right or each side:

$$A^{-1}=E_n\cdots E_2{E}_{1}I$$

This means $A^{-1}$ is obtained by doing the row operations on $I$.

Example

Let's find the inverse of

$$A=\left[\begin{array}{cc}1& 2\\ 3& 4\end{array}\right]$$

A = Matrix([[1,2],[3,4]])
I = eye(2)
aug_mat = A.row_join(I)
aug_mat

$\left[\begin{array}{cccc}1& 2& 1& 0\\ 3& 4& 0& 1\end{array}\right]$

A1 = aug_mat.elementary_row_op("n->n+km",row=1,k=-3,row1=0)
A1

$\left[\begin{array}{cccc}1& 2& 1& 0\\ 0& -2& -3& 1\end{array}\right]$

A2 = A1.elementary_row_op("n->n+km",row=0,k=1,row1=1)
A3 = A2.elementary_row_op("n->kn",row=1,k=-1/2)
A3

$\left[\begin{array}{cccc}1& 0& -2& 1\\ 0& 1& 1.5& -0.5\end{array}\right]$

Therefore, we have

$$A^{-1}=\left[\begin{array}{cc}-2& 1\\ 1.5& -0.5\end{array}\right]$$

We can verify we found the inverse by checking that $A{A}^{-1}=I$.

A_inv = A3[:,2:]
A*A_inv

$\left[\begin{array}{cc}1.0& 0\\ 0& 1.0\end{array}\right]$

Nullspace

The nullspace of a matrix can be computed with .nullspace. This returns an orthogonal basis for the nullspace.

A = Martix([[1,0,0,1],[1,0,1,0]])
A.nullspace()

$[\left[\begin{array}{c}0\\ 1\\ 0\\ 0\end{array}\right],\left[\begin{array}{c}-1\\ 0\\ 1\\ 1\end{array}\right]]$