Solving Linear Systems with Direct Methods#

Teng-Jui Lin

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

Matrix inversion#

Given a system \(\mathbf{A}\mathbf{x} = \mathbf{b}\), we technically have

\[\begin{split} \begin{aligned} \mathbf{A}\mathbf{x} &= \mathbf{b} \\ \mathbf{A}^{-1}\mathbf{A}\mathbf{x} &= \mathbf{A}^{-1}\mathbf{b} \\ \mathbf{x} &= \mathbf{A}^{-1}\mathbf{b} \\ \end{aligned} \end{split}\]

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(\mathbf{A}) = 0\), or \(\mathrm{cond}(\mathbf{A}) \sim 10^{15}\).

Problem Statement. Use matrix inversion to solve the linear system \(\mathbf{Ax=b}\) where

\[\begin{split} \mathbf{A} = \begin{bmatrix} 1 & 0 & 0 \\ 1 & -1 & 0 \\ 2 & -1 & 3 \\ \end{bmatrix}, \mathbf{x} = \begin{bmatrix} x \\ y \\ z \\ \end{bmatrix}, \mathbf{b} = \begin{bmatrix} 5 \\ -2 \\ 9 \\ \end{bmatrix} \end{split}\]
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 \(\mathbf{A}\mathbf{x} = \mathbf{b}\).

We can decompose the matrix \(\mathbf{A}\) into an upper triangular matrix \(\mathbf{L}\) and lower triangular matrix \(\mathbf{U}\):

\[ \mathbf{A} = \mathbf{L}\mathbf{U} \]

so we have

\[\begin{split} \begin{aligned} \mathbf{A}\mathbf{x} &= \mathbf{b} \\ \mathbf{L}\mathbf{U}\mathbf{x} &= \mathbf{b} \\ \mathbf{L}\mathbf{y} &= \mathbf{b} \end{aligned} \end{split}\]

where \(\mathbf{U}\mathbf{x} = \mathbf{y}\). We can first solve for \(\mathbf{y}\), then solve for \(\mathbf{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 \(\mathbf{A}\) is used to solve many different \(\mathbf{b}\).

Problem Statement. Use LU decomposition to solve the linear system \(\mathbf{Ax=b}\) where

\[\begin{split} \mathbf{A} = \begin{bmatrix} 1 & 0 & 0 \\ 1 & -1 & 0 \\ 2 & -1 & 3 \\ \end{bmatrix}, \mathbf{x} = \begin{bmatrix} x \\ y \\ z \\ \end{bmatrix}, \mathbf{b} = \begin{bmatrix} 5 \\ -2 \\ 9 \\ \end{bmatrix} \end{split}\]
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 \(\mathbf{P}\) to change column order, reducing numerical error. In this case, we LU decompose \(\mathbf{PA}\):

\[\begin{split} \begin{aligned} \mathbf{PAx} &= \mathbf{Pb} \\ \mathbf{LUx} &= \mathbf{Pb} \\ \mathbf{L}\mathbf{y} &= \mathbf{Pb} \end{aligned} \end{split}\]

where \(\mathbf{U}\mathbf{x} = \mathbf{y}\). We can then first solve for \(\mathbf{y}\), then solve for \(\mathbf{x}\).

Problem Statement. Use LU decomposition with partial pivoting to solve the linear system \(\mathbf{Ax=b}\) where

\[\begin{split} \mathbf{A} = \begin{bmatrix} 1 & 0 & 0 \\ 1 & -1 & 0 \\ 2 & -1 & 3 \\ \end{bmatrix}, \mathbf{x} = \begin{bmatrix} x \\ y \\ z \\ \end{bmatrix}, \mathbf{b} = \begin{bmatrix} 5 \\ -2 \\ 9 \\ \end{bmatrix} \end{split}\]
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 \(\mathbf{Ax=b}\).

Problem Statement. Use scipy.linalg.solve() to solve the linear system \(\mathbf{Ax=b}\) where

\[\begin{split} \mathbf{A} = \begin{bmatrix} 1 & 0 & 0 \\ 1 & -1 & 0 \\ 2 & -1 & 3 \\ \end{bmatrix}, \mathbf{x} = \begin{bmatrix} x \\ y \\ z \\ \end{bmatrix}, \mathbf{b} = \begin{bmatrix} 5 \\ -2 \\ 9 \\ \end{bmatrix} \end{split}\]
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.])