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

Question

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.

Answer 1

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.

Answer 2

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.

Answer 3

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.

Answer 4

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.