Efficient solver for a symmetric tridiagonal system where the upper/lower diagonals are offset

I'm looking for an efficient way to solve a symmetric tridiagonal system $Mx = d$, where the upper and lower diagonals of $M$ are offset from the main diagonal by $k$ rows/columns:

$$\begin{bmatrix} {b_1 } & {0 } & \ldots & {0 } & {a_1 } & {0 } & \ldots & { 0} \\ {0 } & {b_2 } & {0 } & { } & {0 } & {a_2 } & {0 } & \vdots \\ \vdots & {0 } & {b_3 } & {0 } & { } & {0 } & \ddots & {0 } \\ {0 } & { } & {0 } & {b_4 } & {0 } & { } & { } & {a_{n-k}} \\ {a_1 } & {0 } & { } & {0 } & {b_5 } & {0 } & { } & {0 } \\ {0 } & {a_2 } & {0 } & { } & {0 } & {b_6 } & { } & \vdots \\ \vdots & { } & \ddots & { } & { } & { } & \ddots & {0 } \\ {0 } & \ldots & {0 } & {a_{n-k}} & {0 } & \ldots & {0 } & {b_n } \\ \end{bmatrix} \cdot \begin{bmatrix} {x_ 1 } \\ {x_ 2 } \\ {x_ 3 } \\ {x_ 4 } \\ {x_ 5 } \\ {x_ 6 } \\ \vdots \\ {x_n } \\ \end{bmatrix} = \begin{bmatrix} {d_ 1 } \\ {d_ 2 } \\ {d_ 3 } \\ {d_ 4 } \\ {d_ 5 } \\ {d_ 6 } \\ \vdots \\ {d_n } \\ \end{bmatrix}$$

I'm aware of the TDMA algorithm for solving non-symmetric tridiagonal systems, but I haven't come across an existing implementation that allows the upper and lower diagonals to be arbitrarily offset from the main diagonal.

In addition, I'm wondering whether it would be possible to take advantage of the symmetry of the matrix to solve the system faster than with standard TDMA.

Please excuse the naivety of my question!

Here's a bit more information about $M$:

• $M$ is the Hessian matrix for a large nonlinear system I'm trying to solve via an interior point method
• It is always positive definite
• It is constructed like this: $M = B^{T} \cdot A \cdot B + C$, where $A$ and $C$ are $(n-k, n-k)$ and $(n, n)$ diagonal matrices respectively, and $B$ is an $(n-k, n)$ upper bidiagonal matrix.

I don't think there's much more I can tell you without getting into to much depth, but if you're interested I'm trying to implement the algorithm described in this paper.

• Welcome to SciComp! Does $M$ have any additional structure? (For instance, does it come from a PDE discretization? Where does it come from?) You might be able to use suitable preconditioned iterative methods to solve your problem; these methods would require less memory and would probably scale better. In some cases, destroying the band structure via permutation will yield more favorable matrix orderings for sparse direct solvers. Apr 19, 2014 at 2:10
• @GeoffOxberry I've updated my question with a bit more information about $M$ Apr 19, 2014 at 8:06

When analyzing the matrix structure, one can see that

• column i is connected to column j iff i == j (mod k)
• row i is connected to row j iff i == j (mod k)
• there are at most 2 non-zero elements per row

That leads to the idea of sorting your row and column indices by cosets (w.r.t. "modulo k").

One possible way of sorting by "modulo k" cosets is given by the following permutation myIndex:

n = 7; % example
k = 4;

A = -2*diag(ones(n,1)) + diag(ones(n-k,1),k) + diag(ones(n-k,1),-k);

myIndex = [];
ctr = 1;
subsystemCount = 1;
for p=1:n
myIndex = [myIndex ctr];
ctr = ctr + k;
if ctr > n
subsystemCount = subsystemCount + 1;
if mod(n,k) > 0
ctr = ctr - n;
else
ctr = ctr - n + 1;
end
end
end

B = A(myIndex,myIndex) % tridiagonal system
subsystemCount % number of independent sub-systems

A =
-2     0     0     0     1     0     0
0    -2     0     0     0     1     0
0     0    -2     0     0     0     1
0     0     0    -2     0     0     0
1     0     0     0    -2     0     0
0     1     0     0     0    -2     0
0     0     1     0     0     0    -2

B =
-2     1     0     0     0     0     0
1    -2     0     0     0     0     0
0     0    -2     1     0     0     0
0     0     1    -2     0     0     0
0     0     0     0    -2     1     0
0     0     0     0     1    -2     0
0     0     0     0     0     0    -2


The resulting matrix is tridiagonal and for each coset you get an independent system.

You can then solve all sub systems independently, using an algorithm for solving (symmetric) tridiagonal systems.

Solving the systems independently shouldn’t really be necessary, but it’s interesting to note that you don’t have a single linear system but multiple independent ones for a matrix like that.

After a bit more thought, I think it's actually fairly simple to modify TDMA to deal with arbitrary diagonal offsets:

Forward sweep:

\begin{align} c'_i & = \begin{cases} \cfrac{c_i}{b_i} & ; & i=1,2,\dots,k \\ \cfrac{c_i}{b_i-c'_{i-k} a_i} & ; & i=k+1,k+2,\dots,n-k \end{cases} \\ \\ d'_i &= \begin{cases} \cfrac{d_i}{b_i} & ; & i=1,2,\dots,k \\ \cfrac{d_i-d'_{i-k}a_i}{b_i-c'_{i-k}a_i} & ; & i=k+1,k+2,\dots,n. \\ \end{cases} \end{align}

Back substitution:

\begin{align} x_n &= d'_n\\ x_i &= d'_i-c'_ix_{i+k} & ; {}& i=n-1,n-2,\dots,1. \end{align}

I've written a quick implementation in Cython:

cimport numpy as cnp

# NB: b and d are modified in place
cpdef TDMA_offset(const double[:] a, double[:] b, const double[:] c, double[:] d):
cdef size_t i, k, m, n
m = a.shape[0]
n = b.shape[0]
k = n - m
# Forward sweep
for i in range(m):
d[i+k] -= d[i]*a[i]/b[i]
b[i+k] -= c[i]*a[i]/b[i]
# Back substitution
for i in range(m-1, -1, -1):
d[i] -= d[i+k]*c[i]/b[i+k]
for i in range(n):
d[i] /= b[i]


This approach still doesn't make use of the fact that $M$ is symmetric, and a parallel implementation would also be very nice.

Following the comment by Brian Borchers, I made a simple Python test script:

import numpy as np
from cymodules import _solvers
from scipy import sparse

def test_trisolve(n=10, k=5):
a = np.random.randn(n - k)
b = np.random.randn(n)
c = np.random.randn(n - k)
x = np.random.randn(n)
M = sparse.diags((a, b, c), (-k, 0, k), format='csc')
d = M.dot(x)

# solve for x in place
x_hat = d.copy()
_solvers.TDMA_offset(a, b, c, x_hat)

print "||x - x_hat|| = ", np.linalg.norm(x - x_hat)


I've tested my solver for a variety of values of $n$ and $k$, including cases where $k \geq (n - 1) / 2$:

test_trisolve(n=10, k=3)
# ||x - x_hat|| =  1.67369803253e-15

test_trisolve(n=10, k=6)
# ||x - x_hat|| =  6.85231143806e-16


It seems to be working OK.

• Under the right circumstances (see en.wikipedia.org/wiki/… ) you could specialize the LDL decomposition (see en.wikipedia.org/wiki/Cholesky_decomposition#LDL_decomposition ) the same way you specialized Gaussian elimination for the unsymmetric case. Apr 19, 2014 at 13:09
• I think this breaks if $k \geq (n-1)/2$. Apr 19, 2014 at 15:31
• @BrianBorchers see my update - my Cython implementation seems to work OK in such cases. Perhaps you spotted a mistake in my presentation of the algorithm? Apr 19, 2014 at 18:54
• @JuanM.BelloRivas Sorry, perhaps I wasn't very clear in my description. I'm not actually doing any LDL decomposition in this case - I am directly constructing $M$ from $A$, $B$ and $C$. Apr 19, 2014 at 19:21
• @ali_m you were clear, however TDMA can be obtained by specializing LU decomposition (that is, Gaussian elimination, as in cfm.brown.edu/people/gk/chap6/node13.html ). My suggestion is to follow the same route but starting with LDL instead of LU. Apr 19, 2014 at 20:38

Unfortunately, your matrix isn't tridiagonal in the conventional sense and it can't be turned into a tridiagonal matrix by row/column permutations.

One conventional way of dealing with such a system would be to treat the matrix as a banded matrix and use a factorization routine designed to handle banded matrices. Another approach would be to use a general purpose sparse matrix factorization routine.

One posibility I can think of is by solving many small $2 \times 2$ matrices:

$\begin{bmatrix} {b_1} & {a_1}\\ {a_1} & {b_k}\\ \end{bmatrix} \cdot \begin{bmatrix} {x_1}\\ {x_k}\\ \end{bmatrix} = \begin{bmatrix} {y_1}\\ {y_k}\\ \end{bmatrix}$

Something like,

for all i solve:
$\begin{bmatrix} {b_i} & {a_i}\\ {a_i} & {b_{i+k-1}}\\ \end{bmatrix} \cdot \begin{bmatrix} {x_i}\\ {x_{i+k-1}}\\ \end{bmatrix} = \begin{bmatrix} {y_i}\\ {y_{i+k-1}}\\ \end{bmatrix}$

I think you can make some use of the fact that they are symmetric.

• In the general case, these $2 \times 2$ systems of equations will be linked together. Apr 19, 2014 at 4:07
• Well that's in the general case of a general matrix but unless I misunderstood the structure of the matrix has 2 non zero elements per row and per column. Apr 19, 2014 at 4:10
• I was thinking of the "psychologically tridiagonal" case with wrap around. e.g. if we have a 5 by 5 matrix in which the off diagonal elements occur in the second superdiagonal and the second subdiagonal, then A has the sparsity pattern [X 0 X X 0; 0 X 0 X X; X 0 X 0 X; X X 0 X 0;0 X X 0 X] Looking at the original poster's question that seems not to be the case. Apr 19, 2014 at 4:40
• However, even in the original poster's version, our 5 x 5 example would have the sparsity pattern A=[X 0 X 0 0; 0 X 0 X 0; X 0 X 0 X; 0 X 0 X 0; 0 0 X 0 X]. What would the 2 by 2 systems of equations be in this case? Apr 19, 2014 at 4:41