4.8. Faster Methods for Solving $Ax=b$ for Tridiagonal and Banded matrices, and Strict Diagonal Dominance

Reference:

Section 6.6 Special Types of Matrices in [Burden et al., 2016], the sub-sections on Band Matrices and Tridiagonal Matrices.

Tridiagonal systems

Differential equations often lead to the need to solve systems of equations $Tx=b$ where the matrix $T$ has this speical form:

Definition 4.8 (Tridiagonal matrix)

A matrix $T$ is tridiagonal if the only non-zero elements are on the main diagonal and the diagonals adjacent to it on either side, so that $T_{i,j}=0$ if $|i-j|>1$. That is, the system looks like:

$$\begin{array}{r}Tx=\left[\begin{array}{cc}{d}_1& u_1\\ l_1& d_2& u_2& \\ & l_2& d_3& u_3\\ & & \ddots & \ddots & \ddots \\ & & & l_{n-2}& d_{n-1}& u_{n-1}\\ & & & & l_{n-1}& d_n\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right]=\left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\end{array}\right]\end{array}$$

with all “missing” entries being zeros. The notation used here suggests one efficient way to store such a matrix: as three 1D arrays $d$, $l$ and $u$.

(Such equations also arise in other important situations, such as spline interpolation)

It can be verified that LU factorization preserves all the non-zero values, so that the Doolittle algorithm — if it succeeds without any division by zero — gives $T=LU$ with the form

$$\begin{array}{r}L=\left[\begin{array}{c}1\\ L_1& 1\\ & L_2& 1\\ & & \ddots & \ddots \\ & & & L_{n-2}& 1\\ & & & & L_{n-1}& 1\end{array}\right],U=\left[\begin{array}{cc}{D}_1& u_1\\ & D_2& u_2\\ & & D_3& u_3\\ & & & \ddots & \ddots \\ & & & & D_{n-1}& u_{n-1}\\ & & & & & D_n\end{array}\right]\end{array}$$

Note that the first non-zero element in each column is unchanged, as with a full matrix, but now it means that the upper diagonal elements $u_i$ are unchanged.

Again, one way to describe and store this information is with just the two new 1D arrays $L$ and $D$, along with the unchanged array $u$.

Algorithms

Algorithm 4.11 (LU factorization)

$D_1=d_1$
for i from 2 to n
$L_{i-1}=l_{i-1}/D_{i-1}$
$D_i=d_i-L_{i-1}{u}_{i-1}$
end

Algorithm 4.12 (Forward substitution)

$c_1=b_1$
for i from 2 to n
$c_i=b_i-L_{i-1}{c}_{i-1}$
end

Algorithm 4.13 (Backward substitution)

$x_n=c_n/D_n$
for i from n-1 down to 1
$x_i=(c_i-u_i{x}_{i+1})/D_i$
end

Generalizing to banded matrices

As we have seen, approximating derivatives to higher order of accuracy and approximating derivatives of order greater than two requires more than three nodes, but the locations needed are all close to the ones where the derivative is being approximated. For example, the simplest symmetric approximation of the fourth derivative $D^{4}f(x)$ used values from $f(x-2h)$ to $f(x+2h)$. Then row $i$ of the corresponding matrix has all its non-zero elements at locations $(i,i-2)$ to $(i,i+2)$: the non-zero elements lie in the narrow “band” where $|i-j|\le 2$, and thus on five “central” diagonals.

This is a penta-digonal matrix, and an example of the larger class of banded matrices: ones in which all the non-zero elements have indices $-p\le j-i\le q$ for $p$ and $q$ smaller than $n$ — usually far smaller; $p=q=2$ for a penta-digonal matrix.

Let us recap the general Doolittle algorithm for computing an LU factorization:

Algorithm 4.14 (Doolittle algorithm for computing an LU factorization)

The top row is unchanged:
for j from 1 to n
$$ $u_{1,j}=a_{1,j}$
end
The left column requires no sums:
for i from 2 to n
$$ $l_{i,1}=a_{i,1}/u_{1,1}$
end

The main loop: for k from 2 to n
$$ for j from k to n
$$ $u_{k,j}=a_{k,j}-\sum _{s=1}^{k-1}{l}_{k,s}{u}_{s,j}$
$$ end
$$ for i from k+1 to n
$$ $l_{i,k}=(a_{i,k}-\sum _{s=1}^{k-1}{l}_{i,s}{u}_{s,k})/u_{k,k}$
$$ end
end

Eliminating redundant calculation in the above

With a banded matrix, many of the entries at right are zero, particularly in the two sums, which is where most of the operations are. Thus we can rewrite, exploiting the fact that all elements with indices $j-i<-p$ or $j-i>q$ are zero. To start with, the top diagonal is not modified, as already noted for the tridiagonal case: $u_{k,k+q}=a_{k,k+q}$ for $1\le k\le n-q$.

Algorithm 4.15 (LU factorization of a banded matrix)

The top row is unchanged:
for j from 1 to 1+q
$$ $u_{1,j}=a_{1,j}$
end
The top non-zero diagonal is unchanged:
for k from 1 to n - q
$$ $u_{k,k+q}=a_{k,k+q}$
end
The left column requires no sums:
for i from 2 to 1+p
$$ $l_{i,1}=a_{i,1}/u_{1,1}$
end
The main loop:
for k from 2 to n
$$ for j from k to min(n, k+q-1)
$$ $u_{k,j}=a_{k,j}-{\displaystyle \sum _{s=max(1,k-p,j-q)}^{k-1}{l}_{k,s}{u}_{s,j}}$
$$ end
$$ for i from k+1 to min(n,k+p-1)
$$ $l_{i,k}=(a_{i,k}-{\displaystyle \sum _{s=max(1,i-p,k-q)}^{k-1}{l}_{i,s}{u}_{s,k}})/u_{k,k}$
$$ end
end

It is common for a banded matrix to have equal band-width on either side, $p=q$, as with tridiagonal and pentadiagonal matrices. Then the algorithm is somewhat simpler:

Algorithm 4.16 (LU factorization of a banded matrix, $p=q$)

The top row is unchanged:
for j from 1 to 1+p
$$ $u_{1,j}=a_{1,j}$
end
The top non-zero diagonal is unchanged:
for k from 1 to n - p
$$ $u_{k,k+p}=a_{k,k+p}$
end
The left column requires no sums:
for i from 2 to 1+p
$$ $l_{i,1}=a_{i,1}/u_{1,1}$
end

The main loop:
for k from 2 to n
$$ for j from k to min(n, k+p-1)
$$ $u_{k,j}=a_{k,j}-{\displaystyle \sum _{s=max(1,j-p)}^{k-1}{l}_{k,s}{u}_{s,j}}$
$$ end
$$ for i from k+1 to min(n,k+p)
$$ $l_{i,k}=(a_{i,k}-{\displaystyle \sum _{s=max(1,i-p)}^{k-1}{l}_{i,s}{u}_{s,k}})/u_{k,k}$
$$ end
end

Strict diagonal dominance helps again

These algorithms for banded matrices do no pivoting, and that is highly desirable, because pivoting creates non-zero elements outside the “band” and so can force one back to the general algorithm. Fortunately, we have seen one case where this is fine: the matrix being either row-wise or column-wise strictly diagonally dominant.