Solving Linear Systems with Direct Methods

Solving Linear Systems with Direct Methods#

Teng-Jui Lin

Content adapted from UW AMATH 301, Beginning Scientific Computing, in Spring 2020.

Solving linear systems with direct methods
- Matrix inversion by scipy.linalg.inv()
- LU decomposition with partial pivoting by scipy.linalg.lu()
scipy implementations
- Solving linear systems by scipy.linalg.solve()

Matrix inversion#

Given a system $A x = b$ , we technically have

\begin{array}{r} \begin{aligned} A x & = b \\ A^{- 1} A x & = A^{- 1} b \\ x & = A^{- 1} b \end{aligned} \end{array}

However, we do not use matrix inversion because it is slow, have large error, and not all matrices have an inverse. Singular (non-invertable) matrices have $det (A) = 0$ , or $cond (A) \sim 10^{15}$ .

Problem Statement. Use matrix inversion to solve the linear system $Ax = b$ where

\begin{array}{r} A = [\begin{array}{c} 1 & 0 & 0 \\ 1 & - 1 & 0 \\ 2 & - 1 & 3 \end{array}], x = [\begin{array}{c} x \\ y \\ z \end{array}], b = [\begin{array}{c} 5 \\ - 2 \\ 9 \end{array}] \end{array}

import numpy as np
import scipy
from scipy import linalg

A = np.array([[1, 0, 0], [1, -1, 0], [2, -1, 3]])
b = np.array([5, -2, 9])

# solve by inversion
x = scipy.linalg.inv(A) @ b
x

array([5., 7., 2.])

LU decomposition#

Goal: Solve a system of form $A x = b$ .

We can decompose the matrix $A$ into an upper triangular matrix $L$ and lower triangular matrix $U$ :

A = L U

so we have

\begin{array}{r} \begin{aligned} A x & = b \\ L U x & = b \\ L y & = b \end{aligned} \end{array}

where $U x = y$ . We can first solve for $y$ , then solve for $x$ . Note that upper triangular and lower triangular matrices are easy to solve, so LU decomposition is faster than Gaussian elimination.

LU decomposition is advantageous when the same $A$ is used to solve many different $b$ .

Problem Statement. Use LU decomposition to solve the linear system $Ax = b$ where

\begin{array}{r} A = [\begin{array}{c} 1 & 0 & 0 \\ 1 & - 1 & 0 \\ 2 & - 1 & 3 \end{array}], x = [\begin{array}{c} x \\ y \\ z \end{array}], b = [\begin{array}{c} 5 \\ - 2 \\ 9 \end{array}] \end{array}

A = np.array([[1, 0, 0], [1, -1, 0], [2, -1, 3]])
b = np.array([5, -2, 9])

# solve by LU decomposition
L, U = scipy.linalg.lu(A, permute_l=True)
y = scipy.linalg.solve(L, b)
x = scipy.linalg.solve(U, y)
x

array([5., 7., 2.])

LU decomposition with partial pivoting#

Partial pivoting uses permutation matrix $P$ to change column order, reducing numerical error. In this case, we LU decompose $PA$ :

\begin{array}{r} \begin{aligned} PAx & = Pb \\ LUx & = Pb \\ L y & = Pb \end{aligned} \end{array}

where $U x = y$ . We can then first solve for $y$ , then solve for $x$ .

Problem Statement. Use LU decomposition with partial pivoting to solve the linear system $Ax = b$ where

\begin{array}{r} A = [\begin{array}{c} 1 & 0 & 0 \\ 1 & - 1 & 0 \\ 2 & - 1 & 3 \end{array}], x = [\begin{array}{c} x \\ y \\ z \end{array}], b = [\begin{array}{c} 5 \\ - 2 \\ 9 \end{array}] \end{array}

A = np.array([[1, 0, 0], [1, -1, 0], [2, -1, 3]])
b = np.array([5, -2, 9])

# solve by LU decomposition with partial pivoting
P, L, U = scipy.linalg.lu(A)
y = scipy.linalg.solve(L, P@b)
x = scipy.linalg.solve(U, y)
x

array([5., 7., 2.])

Solving linear systems with `scipy`#

scipy.linalg.solve() directly solves the linear system $Ax = b$ .

Problem Statement. Use scipy.linalg.solve() to solve the linear system $Ax = b$ where

\begin{array}{r} A = [\begin{array}{c} 1 & 0 & 0 \\ 1 & - 1 & 0 \\ 2 & - 1 & 3 \end{array}], x = [\begin{array}{c} x \\ y \\ z \end{array}], b = [\begin{array}{c} 5 \\ - 2 \\ 9 \end{array}] \end{array}

A = np.array([[1, 0, 0], [1, -1, 0], [2, -1, 3]])
b = np.array([5, -2, 9])

# solve by scipy
x = scipy.linalg.solve(A, b)
x

array([5., 7., 2.])