Thomas algorithm for banded matrices

General Usage
1

Dear All,

I am trying to solve linear systems Ax = rhs where A has a tridiagonal shape. When the type of A is tridiagonal, one can use the thomas algo to solve it. I find it faster than defining A as a Tridiagonal and calling \. Indeed, the actual implementation of \ uses a LU decomposition followed by a solve. Maybe the choice of not relying on \ for Tridiagonal is numerical stability of the algorithm.

I want to do something slighly different. I want to solve the same problem for a banded matrix with a diagonal band and two other bands. When the bandwith is 1, one recovers the case above of the Thomas algo.

I have the MWE below, but it does not work. As you can see, I first compute the LU decomposition of the matrix and then solve the linear systems. I have two issues I cannot solve :frowning:

  • Firstly despite my analytical formula, my LU decomposition seems wrong.
  • Secondly, I am able to invert L but I am unable to invert U which I am ashamed of.

If anyone has an idea, Iā€™d be delighted.

Thank you a lot.

using LinearAlgebra, SparseArrays
function thomas2bands!(a::vec, b::vec, c::vec, rhs::vec) where {T, vec <: Vector{T}}
	@assert length(a) == length(c)
	n = length(b)
	p = length(a)
	# we first derive the LU decomposition
	# we have L and U written 8 lines below, L = [l, 1] and U = [v, c]
	v = copy(b)
	l = a ./ v[1:p]
	for ii = n-p+1:n
		v[ii] = b[ii] - l[ii-n+p] * c[ii-n+p]
	end
	# up to here, this has been check, the LU decomposition seems OK
	@show n-p
	L = spdiagm(-n+p => l, 0 => fill(1.0, n))
	U = spdiagm(0 => v, n-p => c)
	A = spdiagm(-n+p => a, 0 => b, n-p => c)
	println("-> A - LU = ", norm(L * U - A, Inf64))

	# we want to solve LUx = rhs
	# we first solve Ly = rhs
	y = similar(rhs)
	for ii=1:n-p
		y[ii] = rhs[ii]
	end
	for ii=n-p+1:n
		y[ii] = rhs[ii] - l[ii-n+p] * y[ii-n+p]
	end
	println("-> L \\ rhs = ",norm(y -  (L \ rhs), Inf64))

	# we solve Ux = y
	x = similar(rhs)
	for ii = n-p+1:n
		x[ii] = y[ii] / v[ii]
	end
	for ii=1:n-p
		x[ii] = y[ii] / v[ii] - x[ii+n-p] * c[ii] / v[ii]
	end
	println("-> y - U * x   = ",norm((y -  (U * x)), Inf64))
	println("-> U \\ y   = ",norm(x -  (U \ y),Inf64))
	return x
end
	sol = @time thomas2bands!(rand(10^5-10^3), 2 .+ rand(10^5), rand(10^5-10^3), rand(10^5) )
1 Like
4

I think the problem is that in the analytical case you consider to build the full formula, you assume that there is no column where both off diagonal terms are non-zero. For instance consider the 3d column in

5

Thank you for your reply! When I multiply your matrices, I get

 1  0  1  0  0
 0  1  0  1  0
 1  0  2  0  1
 0  1  0  2  0
 0  0  1  0  2

It seems it has the correct sparsity pattern. Did you look at my Maple document?

6

Ha sorry I was checking what I wrote on a piece of paper as you typed your response and realised thatā€™s not it (sorry) and yes I did look at the maple document I thought it was missing a case and that it could explain the problem.

Edit: there still seems to be a problem with the lower diagonal. When doing

sol, R = @time thomas2bands!(rand(3), 2 .+ rand(5), rand(3), rand(5) )

where R is Matrix(A-L*U) I get

5Ɨ5 Array{Float64,2}:
 0.0  0.0  0.0        0.0  0.0
 0.0  0.0  0.0        0.0  0.0
 0.0  0.0  0.0        0.0  0.0
 0.0  0.0  0.0        0.0  0.0
 0.0  0.0  0.0190783  0.0  0.0

however when going with 6, 3 then itā€™s fine. with 7 and 5 the off diagonal is wrong for the last 3 elements.

7

I agree!! Thatā€™s my first mistake I canā€™t solve. This is strange because the v part seems fine despite involving l which seems to be wrongā€¦

You can check this by doing lu(A), the sparsity is the same as I assumed for L and Uā€¦

BTW: I updated the document with the formula.

8

Actualy thatā€™s just what I did and it doesnā€™t seem to be the case? (this is for 7/5)

julia> L
6Ɨ6 Array{Float64,2}:
 1.0       0.0       0.0  0.0  0.0  0.0
 0.0       1.0       0.0  0.0  0.0  0.0
 0.0       0.0       1.0  0.0  0.0  0.0
 0.0       0.0       0.0  1.0  0.0  0.0
 0.823648  0.0       0.0  0.0  1.0  0.0
 0.0       0.910357  0.0  0.0  0.0  1.0

julia> luL
7Ɨ7 Array{Float64,2}:
 1.0       0.0       0.0        0.0       0.0      0.0  0.0
 0.0       1.0       0.0        0.0       0.0      0.0  0.0
 0.220607  0.0       1.0        0.0       0.0      0.0  0.0
 0.0       0.209361  0.0        1.0       0.0      0.0  0.0
 0.0       0.0       0.0931903  0.0       1.0      0.0  0.0
 0.0       0.0       0.0        0.324195  0.0      1.0  0.0
 0.0       0.0       0.0        0.0       0.06044  0.0  1.0

julia> U
6Ɨ6 Array{Float64,2}:
 0.164566  0.0       0.0      0.0       0.361828   0.0      
 0.0       0.177329  0.0      0.0       0.0        0.973216 
 0.0       0.0       0.27888  0.0       0.0        0.0      
 0.0       0.0       0.0      0.203477  0.0        0.0      
 0.0       0.0       0.0      0.0       0.0423017  0.0      
 0.0       0.0       0.0      0.0       0.0        0.0682693

julia> luU
7Ɨ7 Array{Float64,2}:
 2.65545  0.0      0.353129  0.0       0.0       0.0       0.0      
 0.0      2.57589  0.0       0.767602  0.0       0.0       0.0      
 0.0      0.0      2.79038   0.0       0.043141  0.0       0.0      
 0.0      0.0      0.0       2.80709   0.0       0.267985  0.0      
 0.0      0.0      0.0       0.0       2.76367   0.0       0.0668464
 0.0      0.0      0.0       0.0       0.0       2.38242   0.0      
 0.0      0.0      0.0       0.0       0.0       0.0       2.05833

Edit: hmmm the dimensions are not quite rightā€¦ could that be it?

9

lu(Matrix(A)) I guess

Also, there is a mistake in the definition of L which I corrected

11

It seems it works when the bandwith is even actually which is quite puzzling. Maybe a mistake with the indices

12

Ok I spent quite some time on this because I thought it was fun. Hereā€™s an algo that actually works, note that the decomposition is actually a bit trickier.

Itā€™s just a bit annoying with the indices floating around


function decomp(A, p)
    n = size(A, 1)

    Ī» = zeros(p)
    Ī¼ = zeros(p)
    v = zeros(n)

    for i=1:n-p
        v[i] = A[i,i]
    end

    for j = n-p+1:n
        Ī¼[j-n+p] = A[j-n+p, j]
    end

    for k=1:Int(floor(n/(n-p)))
        for i=(k*(n-p)+1):min((k+1)*(n-p), n)
            Ī»[i-n+p] = A[i,i-n+p] / v[i-n+p]
            v[i] = A[i, i] - Ī»[i-n+p] * Ī¼[i-n+p]
        end
    end

    v, Ī», Ī¼
end

buildL(Ī», n) = diagm(-n+length(Ī») => Ī») + I

buildU(v, Ī¼, n) = diagm(0=>v, n-length(Ī¼)=>Ī¼)


n, p = 7, 5

v = randn(n)
Ī» = randn(p)
Ī¼ = randn(p)

L = buildL(Ī», n)
U = buildU(v, Ī¼, n)

A = L*U

rv, rĪ», rĪ¼ = decomp(A, p)

norm(rv - v)
norm(rĪ» - Ī»)
norm(rĪ¼ - Ī¼)

Edit: to get to that, I wrote down the equation for rows of L and columns of U and computed what their dot product looks like in general. This gave me:

function rowL(i, Ī», n)
    p = length(Ī»)
    ri = zeros(n)
    for j = 1:n
        ri[j] = (j==i ? 1.0 : 0.0) +
                   ((j == i-n+p) ? Ī»[j] : 0.0)
    end
    ri
end

function colU(j, v, Ī¼)
    n  = length(v)
    p  = length(Ī¼)
    cj = zeros(n)
    for i = 1:n
        cj[i] = (i==j ? v[i] : 0.0) +
                    ((i == j-n+p) ? Ī¼[i] : 0.0)
    end
    cj
end

function LU(i, j, v, Ī», Ī¼)
    p = length(Ī»)
    n = length(v)

    ri = rowL(i, Ī», n)
    cj = colU(j, v, Ī¼)

    gt = dot(ri, cj)

    hand = (i == j ? v[i] + (i ā‰„ n-p+1 ? Ī»[i-n+p]*Ī¼[i-n+p] : 0.0) : 0.0) +
           (j == i-n+p ? Ī»[i-n+p]*v[i-n+p] : 0.0) +
           (i == j-n+p ? Ī¼[i] : 0.0)

#    @show  gt-hand
    return hand
end

The expressions here are not simplified to clarify how I got an expression for a[i, j]. Once I had the LU bit, I reversed the steps and realised that there was a trick where you can only get v and lambda in blocks (as per previous algo)

Edit2 what I mean by ā€œblocksā€ is that the equations allow you to get things from n-p+1 to 2n-2p and then from 2n-2p+1 to 3n-3p etc until you bottom out. But past n-p+1 you cannot in general get all remaining v and lambda in one shot.

I think the Maple document you had doesnā€™t show this because p is too small but if you take the n=7 and p=5 case, you can see how some of this plays out.

Edit3: it may be that the decomposition I suggest can be further simplified, I didnā€™t really try, I just know that itā€™s correct as I tried a big bunch of n and p and the relative error is always around 1e-15.

13

Hereā€™s your code adjusted


using LinearAlgebra, SparseArrays
function thomas2bands!(a::vec, b::vec, c::vec, rhs::vec) where {T, vec <: Vector{T}}
	@assert length(a) == length(c)
	n = length(b)
	p = length(a)
	# we first derive the LU decomposition
	# we have L and U written 8 lines below, L = [l, 1] and U = [v, c]

	Ī» = zeros(p)
    Ī¼ = zeros(p)
    v = zeros(n)

    for i=1:n-p
        v[i] = b[i]
    end

    for k=1:Int(floor(n/(n-p)))
        for i=(k*(n-p)+1):min((k+1)*(n-p), n)
            Ī»[i-n+p] = a[i-n+p] / v[i-n+p]
            v[i] = b[i] - Ī»[i-n+p] * c[i-n+p]
        end
    end

	# up to here, this has been check, the LU decomposition seems OK
	@show n-p
	L = spdiagm(-n+p => Ī», 0 => fill(1.0, n))
	U = spdiagm(0 => v, n-p => c)
	A = spdiagm(-n+p => a, 0 => b, n-p => c)
	println("-> A - LU = ", norm(L * U - A, Inf64))

	# we want to solve LUx = rhs
	# we first solve Ly = rhs
	y = similar(rhs)
	for ii=1:n-p
		y[ii] = rhs[ii]
	end
	for ii=n-p+1:n
		y[ii] = rhs[ii] - Ī»[ii-n+p] * y[ii-n+p]
	end
	println("-> L \\ rhs = ",norm(y -  (L \ rhs), Inf64))

	# we solve Ux = y
	x = similar(rhs)
	for ii = n-p+1:n
		x[ii] = y[ii] / v[ii]
	end
	for ii=1:n-p
		x[ii] = y[ii] / v[ii] - x[ii+n-p] * c[ii] / v[ii]
	end
	println("-> y - U * x   = ",norm((y -  (U * x)), Inf64))
	println("-> U \\ y   = ",norm(x -  (U \ y),Inf64))
	return x
end
	sol = @time thomas2bands!(rand(10^5-10^3), 2 .+ rand(10^5), rand(10^5-10^3), rand(10^5))

This gives

n - p = 1000
-> A - LU = 4.440892098500626e-16

(I didnā€™t check the rest)

14

Thank you a lot for your intense help!

I will digest this.

1 Like
15

Hey

Sorry for coming back late. The following work for me

function thomas2bands!(a, b, c, d; debug = false)
	@assert length(a) == length(c)
	@assert length(b) == length(d)
	debug && println("#"^20)
	n = length(b)
	p = length(a)
	# we first derive the LU decomposition
	# we have L and U written 8 lines below, L = [l, 1] and U = [v, c]
	v = similar(b)
	l = similar(a)
	for i=1:n-p
		v[i] = b[i]
	end

	for i=n-p+1:n
		l[i-n+p] = a[i-n+p] / v[i-n+p]
		v[i] = b[i] - l[i-n+p] * c[i-n+p]
	end
	# up to here, this has been check, the LU decomposition seems OK

	debug && (L = spdiagm(-n+p => l, 0 => fill(1.0, n)))
	debug && (U = spdiagm(0 => v, n-p => c))
	debug && (A = spdiagm(-n+p => a, 0 => b, n-p => c))
	debug && println("-> A - LU = ", norm(L * U - A, Inf64))
end
16

Hello Romain,

I am solving a Richards equation (hydrology) by using a tridiagonal nonlinear set of equations in hydrology by using the Newton-Raphson method. I understand that your code will do the job. Is your posted code the latest version?

Thanks for any help you may provide?

Banded matrice
Banded matrice1608Ɨ501 76.7 KB

17

Hi,

You may be mistaken. I was interested in a generalisation of the Thomas algorithm when there are 5 bands. You can invert your matrix by making it a Tridiagonal and calling \.

1 Like