381 - A20 - Banded Systems

145 days ago by Professor381

John Travis

Mississippi College

When solving many practical problems involving linear systems, it is the case that the resulting matrix A may contain a significant number of zero entries.  This happens often when solving (for example) differential equation problems and when determining the coefficients for cubic splines. These zeros can be exploited to develop algorithms for solving $Ax=b$ using only $O(n)$ calculations and signficantly less storage requirements. These algorithms build off of the matrix factorization techniques that put $A = LU$ and then storing all of the matrices not by horizontal rows and vertical columns but by focusing on "diagonals".

Below, let's consider a case when A is tridiagonal...that is only has non-zero entries on the main diagonal and on the diagonals one above and one below.
n = 7 A = zeros(n,n); for k=1:n A(k,k) = floor(10*rand+1)+10; # let's make A diagonally dominant b(k) = floor(10*rand+1); # we will use this later to finally use the factorization to solve Ax=b end for k = 1:n-1 A(k,k+1) = floor(10*rand+1); A(k+1,k) = floor(10*rand+1); end format short A b = b' 
        
n = 7
A =

   17    8    0    0    0    0    0
    2   17    9    0    0    0    0
    0    1   11    9    0    0    0
    0    0    6   11    8    0    0
    0    0    0    4   13    4    0
    0    0    0    0    1   19    9
    0    0    0    0    0    5   19

b =

    4
    4
    5
   10
    7
    2
    3
n = 7
A =

   17    8    0    0    0    0    0
    2   17    9    0    0    0    0
    0    1   11    9    0    0    0
    0    0    6   11    8    0    0
    0    0    0    4   13    4    0
    0    0    0    0    1   19    9
    0    0    0    0    0    5   19

b =

    4
    4
    5
   10
    7
    2
    3

Notice all of the zeros off of these three diagonals.  Since they are zeros that won't change using Gaussian Elimination then let's not even store them.  To do so, let's just store the values in the three diagonal into three long vectors.  In essence, this means the row values in each of these vectors will be correct but the column points will be messed up.  Indeed, let's store

$A_{k+1,k}$ into band$(k+1,1)$

$A_{k,k}$ into band$(k,2)$

$A_{k,k+1}$ into band$(k,3)$

which would give the following storage for the matrix above:
for k=1:n band(k,2) = A(k,k); end band(1,1)=0; band(n,3)=0; for k = 1:n-1 band(k,3) = A(k,k+1); band(k+1,1) = A(k+1,k); end disp('Here is the original A:') format short A disp('Here it is using banded storage:') band 
        
Here is the original A:
A =

   17    8    0    0    0    0    0
    2   17    9    0    0    0    0
    0    1   11    9    0    0    0
    0    0    6   11    8    0    0
    0    0    0    4   13    4    0
    0    0    0    0    1   19    9
    0    0    0    0    0    5   19

Here it is using banded storage:
band =

    0   17    8
    2   17    9
    1   11    9
    6   11    8
    4   13    4
    1   19    9
    5   19    0
Here is the original A:
A =

   17    8    0    0    0    0    0
    2   17    9    0    0    0    0
    0    1   11    9    0    0    0
    0    0    6   11    8    0    0
    0    0    0    4   13    4    0
    0    0    0    0    1   19    9
    0    0    0    0    0    5   19

Here it is using banded storage:
band =

    0   17    8
    2   17    9
    1   11    9
    6   11    8
    4   13    4
    1   19    9
    5   19    0

Notice that the required memory storage for the matrix A is now $O(3n)$ rather than $O(n^2)$.

We can now modify our existing factorization, forward solving, and backward solving routines to utilize this triadiagonal structure and storage scheme.
# For a tridiagonal system, factor into A=LU with U having 1's on the diagonal # Storage will be mat = [ lower diagonal of L , main diagonal of L, upper diagonal of U] function ans = factorLU_tri(band) n = size(band)(1); band(1,3) = band(1,3)/band(1,2); for step = 2:n-1 band(step,2) = band(step,2)-band(step,1)*band(step-1,3); band(step,3) = band(step,3)/band(step,2); end band(n,2) = band(n,2)-band(n,1)*band(n-1,3); ans = band; endfunction function ans = forwardsolve_tri(band,rhs) # with mat in 3-column storage n = size(rhs)(1); y(1) = rhs(1)/band(1,2); # the first equation for k = 2:n # each successive equation y(k) = (rhs(k)-band(k,1)*y(k-1))/band(k,2); end ans = y'; endfunction function ans = backwardsolve_tri(band,rhs) #with mat in 3-column storage n = size(rhs)(1); x = zeros(n,1); x(n) = rhs(n); # the last equation for k = n-1:-1:1 # each successive equation x(k) = rhs(k)-band(k,3)*x(k+1); # Note that U has one's on the diagonal end ans = x'; endfunction 
       
clear x mat_factored = factorLU_tri(band) disp('Solving using') y = forwardsolve_tri(mat_factored,b) disp('yields') x = backwardsolve_tri(mat_factored,y); x = x' disp('with b-Ax = ') b-A*x 
        
mat_factored =

    0.00000   17.00000    0.47059
    2.00000   16.05882    0.56044
    1.00000   10.43956    0.86211
    6.00000    5.82737    1.37283
    4.00000    7.50867    0.53272
    1.00000   18.46728    0.48735
    5.00000   16.56326    0.00000

Solving using
y =

   0.235294
   0.219780
   0.457895
   1.244581
   0.269246
   0.093720
   0.152832

yields
x =

   0.050496
   0.392695
  -0.308535
   0.889021
   0.258997
   0.019237
   0.152832

with b-Ax = 
ans =

   0
   0
   0
   0
   0
   0
   0
mat_factored =

    0.00000   17.00000    0.47059
    2.00000   16.05882    0.56044
    1.00000   10.43956    0.86211
    6.00000    5.82737    1.37283
    4.00000    7.50867    0.53272
    1.00000   18.46728    0.48735
    5.00000   16.56326    0.00000

Solving using
y =

   0.235294
   0.219780
   0.457895
   1.244581
   0.269246
   0.093720
   0.152832

yields
x =

   0.050496
   0.392695
  -0.308535
   0.889021
   0.258997
   0.019237
   0.152832

with b-Ax = 
ans =

   0
   0
   0
   0
   0
   0
   0

Now, modify the routines above to handle a 5-diagonal system with two above and two below the main diagonal.