Solving Linear Systems

import numpy as np
import matplotlib.pyplot as plt
import scipy.linalg as la
%matplotlib inline

Linear Systems

A linear system of equations is a collection of linear equations

\begin{align} a_{0,0}x_0 + a_{0,1}x_2 + \cdots + a_{0,n}x_n & = b_0 \\ a_{1,0}x_0 + a_{1,1}x_2 + \cdots + a_{1,n}x_n & = b_1 \\ & \vdots \\ a_{m,0}x_0 + a_{m,1}x_2 + \cdots + a_{m,n}x_n & = b_m \\ \end{align}

In matrix notation, a linear system is $A \mathbf{x}= \mathbf{b}$ where

$$ A = \begin{bmatrix} a_{0,0} & a_{0,1} & \cdots & a_{0,n} \\ a_{1,0} & a_{1,1} & \cdots & a_{1,n} \\ \vdots & & & \vdots \\ a_{m,0} & a_{m,1} & \cdots & a_{m,n} \\ \end{bmatrix} \ \ , \ \ \mathbf{x} = \begin{bmatrix} x_0 \\ x_1 \\ \vdots \\ x_n \end{bmatrix} \ \ , \ \ \mathbf{b} = \begin{bmatrix} b_0 \\ b_1 \\ \vdots \\ b_m \end{bmatrix} $$

Gaussian elimination

The general procedure to solve a linear system of equation is called Gaussian elimination. The idea is to perform elementary row operations to reduce the system to its row echelon form and then solve.

Elementary Row Operations

Elementary row operations include:

  1. Add $k$ times row $j$ to row $i$.
  2. Multiply row $i$ by scalar $k$.
  3. Switch rows $i$ and $j$.

Each of the elementary row operations is the result of matrix multiplication by an elementary matrix (on the left).

To add $k$ times row $i$ to row $j$ in a matrix $A$, we multiply $A$ by the matrix $E$ where $E$ is equal to the identity matrix except $E_{i,j} = k$. For example, if $A$ is 3 by 3 and we want to add 3 times row 2 to row 0 (using 0 indexing) then

$$ E = \begin{bmatrix} 1 & 0 & 3 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} $$

Let's verify the calculation:

A = np.array([[1,1,2],[-1,3,1],[0,5,2]])
print(A)
[[ 1  1  2]
 [-1  3  1]
 [ 0  5  2]]
E1 = np.array([[1,0,3],[0,1,0],[0,0,1]])
print(E1)
[[1 0 3]
 [0 1 0]
 [0 0 1]]
E1 @ A
array([[ 1, 16,  8],
       [-1,  3,  1],
       [ 0,  5,  2]])

To multiply $k$ times row $i$ in a matrix $A$, we multiply $A$ by the matrix $E$ where $E$ is equal to the identity matrix except $E_{i,i} = k$. For example, if $A$ is 3 by 3 and we want to multiply row 1 by -2 then

$$ E = \begin{bmatrix} 1 & 0 & 0 \\ 0 & -2 & 0 \\ 0 & 0 & 1 \end{bmatrix} $$

Let's verify the calculation:

E2 = np.array([[1,0,0],[0,-2,0],[0,0,1]])
print(E2)
[[ 1  0  0]
 [ 0 -2  0]
 [ 0  0  1]]
E2 @ A
array([[ 1,  1,  2],
       [ 2, -6, -2],
       [ 0,  5,  2]])

Finally, to switch row $i$ and row $j$ in a matrix $A$, we multiply $A$ by the matrix $E$ where $E$ is equal to the identity matrix except $E_{i,i} = 0$, $E_{j,j} = 0$, $E_{i,j} = 1$ and $E_{j,i} = 1$. For example, if $A$ is 3 by 3 and we want to switch row 1 and row 2 then

$$ E = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \end{bmatrix} $$

Let's verify the calculation:

E3 = np.array([[1,0,0],[0,0,1],[0,1,0]])
print(E3)
[[1 0 0]
 [0 0 1]
 [0 1 0]]
E3 @ A
array([[ 1,  1,  2],
       [ 0,  5,  2],
       [-1,  3,  1]])

Implementation

Let's write function to implement the elementary row operations. First of all, let's write a function called add_rows which takes input parameters $A$, $k$, $i$ and $j$ and returns the NumPy array resulting from adding $k$ times row $j$ to row $i$ in the matrix $A$. If $i=j$, then let's say that the function scales row $i$ by $k+1$ since this would be the result of $k$ times row $i$ added to row $i$.

def add_row(A,k,i,j):
    "Add k times row j to row i in matrix A."
    n = A.shape[0]
    E = np.eye(n)
    if i == j:
        E[i,i] = k + 1
    else:
        E[i,j] = k
    return E @ A

Let's test our function:

M = np.array([[1,1],[3,2]])
print(M)
[[1 1]
 [3 2]]
add_row(M,2,0,1)
array([[7., 5.],
       [3., 2.]])
add_row(M,3,1,1)
array([[ 1.,  1.],
       [12.,  8.]])

Let's write a function called scale_row which takes 3 input parameters $A$, $k$, and $i$ and returns the matrix that results from $k$ times row $i$ in the matrix $A$.

def scale_row(A,k,i):
    "Multiply row i by k in matrix A."
    n = A.shape[0]
    E = np.eye(n)
    E[i,i] = k
    return E @ A
M = np.array([[3,1],[-2,7]])
print(M)
[[ 3  1]
 [-2  7]]
scale_row(M,3,1)
array([[ 3.,  1.],
       [-6., 21.]])
A = np.array([[1,1,1],[1,-1,0]])
print(A)
[[ 1  1  1]
 [ 1 -1  0]]
scale_row(A,5,1)
array([[ 1.,  1.,  1.],
       [ 5., -5.,  0.]])

Let's write a function called switch_rows which takes 3 input parameters $A$, $i$ and $j$ and returns the matrix that results from switching rows $i$ and $j$ in the matrix $A$.

def switch_rows(A,i,j):
    "Switch rows i and j in matrix A."
    n = A.shape[0]
    E = np.eye(n)
    E[i,i] = 0
    E[j,j] = 0
    E[i,j] = 1
    E[j,i] = 1
    return E @ A
A = np.array([[1,1,1],[1,-1,0]])
print(A)
[[ 1  1  1]
 [ 1 -1  0]]
switch_rows(A,0,1)
array([[ 1., -1.,  0.],
       [ 1.,  1.,  1.]])

Examples

Find the Inverse

Let's apply our functions to the augmented matrix $[M \ | \ I]$ to find the inverse of the matrix $M$:

M = np.array([[5,4,2],[-1,2,1],[1,1,1]])
print(M)
[[ 5  4  2]
 [-1  2  1]
 [ 1  1  1]]
A = np.hstack([M,np.eye(3)])
print(A)
[[ 5.  4.  2.  1.  0.  0.]
 [-1.  2.  1.  0.  1.  0.]
 [ 1.  1.  1.  0.  0.  1.]]
A1 = switch_rows(A,0,2)
print(A1)
[[ 1.  1.  1.  0.  0.  1.]
 [-1.  2.  1.  0.  1.  0.]
 [ 5.  4.  2.  1.  0.  0.]]
A2 = add_row(A1,1,1,0)
print(A2)
[[1. 1. 1. 0. 0. 1.]
 [0. 3. 2. 0. 1. 1.]
 [5. 4. 2. 1. 0. 0.]]
A3 = add_row(A2,-5,2,0)
print(A3)
[[ 1.  1.  1.  0.  0.  1.]
 [ 0.  3.  2.  0.  1.  1.]
 [ 0. -1. -3.  1.  0. -5.]]
A4 = switch_rows(A3,1,2)
print(A4)
[[ 1.  1.  1.  0.  0.  1.]
 [ 0. -1. -3.  1.  0. -5.]
 [ 0.  3.  2.  0.  1.  1.]]
A5 = scale_row(A4,-1,1)
print(A5)
[[ 1.  1.  1.  0.  0.  1.]
 [ 0.  1.  3. -1.  0.  5.]
 [ 0.  3.  2.  0.  1.  1.]]
A6 = add_row(A5,-3,2,1)
print(A6)
[[  1.   1.   1.   0.   0.   1.]
 [  0.   1.   3.  -1.   0.   5.]
 [  0.   0.  -7.   3.   1. -14.]]
A7 = scale_row(A6,-1/7,2)
print(A7)
[[ 1.          1.          1.          0.          0.          1.        ]
 [ 0.          1.          3.         -1.          0.          5.        ]
 [ 0.          0.          1.         -0.42857143 -0.14285714  2.        ]]
A8 = add_row(A7,-3,1,2)
print(A8)
[[ 1.          1.          1.          0.          0.          1.        ]
 [ 0.          1.          0.          0.28571429  0.42857143 -1.        ]
 [ 0.          0.          1.         -0.42857143 -0.14285714  2.        ]]
A9 = add_row(A8,-1,0,2)
print(A9)
[[ 1.          1.          0.          0.42857143  0.14285714 -1.        ]
 [ 0.          1.          0.          0.28571429  0.42857143 -1.        ]
 [ 0.          0.          1.         -0.42857143 -0.14285714  2.        ]]
A10 = add_row(A9,-1,0,1)
print(A10)
[[ 1.          0.          0.          0.14285714 -0.28571429  0.        ]
 [ 0.          1.          0.          0.28571429  0.42857143 -1.        ]
 [ 0.          0.          1.         -0.42857143 -0.14285714  2.        ]]

Let's verify that we found the inverse $M^{-1}$ correctly:

Minv = A10[:,3:]
print(Minv)
[[ 0.14285714 -0.28571429  0.        ]
 [ 0.28571429  0.42857143 -1.        ]
 [-0.42857143 -0.14285714  2.        ]]
result = Minv @ M
print(result)
[[ 1.00000000e+00  2.22044605e-16  1.11022302e-16]
 [-4.44089210e-16  1.00000000e+00 -2.22044605e-16]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]

Success! We can see the result more clearly if we round to 15 decimal places:

np.round(result,15)
array([[ 1.,  0.,  0.],
       [-0.,  1., -0.],
       [ 0.,  0.,  1.]])

Solve a System

Let's use our functions to perform Gaussian elimination and solve a linear system of equations $A \mathbf{x} = \mathbf{b}$.

A = np.array([[6,15,1],[8,7,12],[2,7,8]])
print(A)
[[ 6 15  1]
 [ 8  7 12]
 [ 2  7  8]]
b = np.array([[2],[14],[10]])
print(b)
[[ 2]
 [14]
 [10]]

Form the augemented matrix $M$:

M = np.hstack([A,b])
print(M)
[[ 6 15  1  2]
 [ 8  7 12 14]
 [ 2  7  8 10]]

Perform row operations:

M1 = scale_row(M,1/6,0)
print(M1)
[[ 1.          2.5         0.16666667  0.33333333]
 [ 8.          7.         12.         14.        ]
 [ 2.          7.          8.         10.        ]]
M2 = add_row(M1,-8,1,0)
print(M2)
[[  1.           2.5          0.16666667   0.33333333]
 [  0.         -13.          10.66666667  11.33333333]
 [  2.           7.           8.          10.        ]]
M3 = add_row(M2,-2,2,0)
print(M3)
[[  1.           2.5          0.16666667   0.33333333]
 [  0.         -13.          10.66666667  11.33333333]
 [  0.           2.           7.66666667   9.33333333]]
M4 = scale_row(M3,-1/13,1)
print(M4)
[[ 1.          2.5         0.16666667  0.33333333]
 [ 0.          1.         -0.82051282 -0.87179487]
 [ 0.          2.          7.66666667  9.33333333]]
M5 = add_row(M4,-2,2,1)
print(M5)
[[ 1.          2.5         0.16666667  0.33333333]
 [ 0.          1.         -0.82051282 -0.87179487]
 [ 0.          0.          9.30769231 11.07692308]]
M6 = scale_row(M5,1/M5[2,2],2)
print(M6)
[[ 1.          2.5         0.16666667  0.33333333]
 [ 0.          1.         -0.82051282 -0.87179487]
 [ 0.          0.          1.          1.19008264]]
M7 = add_row(M6,-M6[1,2],1,2)
print(M7)
[[1.         2.5        0.16666667 0.33333333]
 [0.         1.         0.         0.1046832 ]
 [0.         0.         1.         1.19008264]]
M8 = add_row(M7,-M7[0,2],0,2)
print(M8)
[[1.         2.5        0.         0.13498623]
 [0.         1.         0.         0.1046832 ]
 [0.         0.         1.         1.19008264]]
M9 = add_row(M8,-M8[0,1],0,1)
print(M9)
[[ 1.          0.          0.         -0.12672176]
 [ 0.          1.          0.          0.1046832 ]
 [ 0.          0.          1.          1.19008264]]

Success! The solution of $Ax=b$ is

x = M9[:,3].reshape(3,1)
print(x)
[[-0.12672176]
 [ 0.1046832 ]
 [ 1.19008264]]

Or, we can do it the easy way...

x = la.solve(A,b)
print(x)
[[-0.12672176]
 [ 0.1046832 ]
 [ 1.19008264]]

scipy.linalg.solve

We are mostly interested in linear systems $A \mathbf{x} = \mathbf{b}$ where there is a unique solution $\mathbf{x}$. This is the case when $A$ is a square matrix ($m=n$) and $\mathrm{det}(A) \not= 0$. To solve such a system, we can use the function scipy.linalg.solve.

The function returns a solution of the system of equations $A \mathbf{x} = \mathbf{b}$. For example:

A = np.array([[1,1],[1,-1]])
print(A)
[[ 1  1]
 [ 1 -1]]
b1 = np.array([2,0])
print(b1)
[2 0]

And solve:

x1 = la.solve(A,b1)
print(x1)
[1. 1.]

Note that the output $\mathbf{x}$ is returned as a 1D NumPy array when the vector $\mathbf{b}$ (the right hand side) is entered as a 1D NumPy array. If we input $\mathbf{b}$ as a 2D NumPy array, then the output is a 2D NumPy array. For example:

A = np.array([[1,1],[1,-1]])
b2 = np.array([2,0]).reshape(2,1)
x2 = la.solve(A,b2)
print(x2)
[[1.]
 [1.]]

Finally, if the right hand side $\mathbf{b}$ is a matrix, then the output is a matrix of the same size. It is the solution of $A \mathbf{x} = \mathbf{b}$ when $\mathbf{b}$ is a matrix. For example:

A = np.array([[1,1],[1,-1]])
b3 = np.array([[2,2],[0,1]])
x3 = la.solve(A,b3)
print(x3)
[[1.  1.5]
 [1.  0.5]]

Simple Example

Let's compute the solution of the system of equations

\begin{align} 2x + y &= 1 \\ x + y &= 1 \end{align}

Create the matrix of coefficients:

A = np.array([[2,1],[1,1]])
print(A)
[[2 1]
 [1 1]]

And the vector $\mathbf{b}$:

b = np.array([1,-1]).reshape(2,1)
print(b)
[[ 1]
 [-1]]

And solve:

x = la.solve(A,b)
print(x)
[[ 2.]
 [-3.]]

We can verify the solution by computing the inverse of $A$:

Ainv = la.inv(A)
print(Ainv)
[[ 1. -1.]
 [-1.  2.]]

And multiply $A^{-1} \mathbf{b}$ to solve for $\mathbf{x}$:

x = Ainv @ b
print(x)
[[ 2.]
 [-3.]]

We get the same result. Success!

Inverse or Solve

It's a bad idea to use the inverse $A^{-1}$ to solve $A \mathbf{x} = \mathbf{b}$ if $A$ is large. It's too computationally expensive. Let's create a large random matrix $A$ and vector $\mathbf{b}$ and compute the solution $\mathbf{x}$ in 2 ways:

N = 1000
A = np.random.rand(N,N)
b = np.random.rand(N,1)

Check the first entries $A$:

A[:3,:3]
array([[0.86338109, 0.60332412, 0.91735417],
       [0.56476439, 0.90068576, 0.85848553],
       [0.35391662, 0.79723924, 0.61917049]])

And for $\mathbf{b}$:

b[:4,:]
array([[0.04358765],
       [0.8484762 ],
       [0.15199503],
       [0.99150553]])

Now we compare the speed of scipy.linalg.solve with scipy.linalg.inv:

%%timeit
x = la.solve(A,b)
30.2 ms ± 511 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
%%timeit
x = la.inv(A) @ b
54.2 ms ± 1.36 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

Solving with scipy.linalg.solve is about twice as fast!

Exercises

Under construction