Cyclic Tridiagonal Matrices

In a recent post, I talked about how the Thomas algorithm can be used to solve a tridiagonal system of equations very efficiently. Recall, that a tridiagonal matrix looks like, \[A = \begin{bmatrix} b_1 & c_1 & 0 & ... & 0\\ a_2 & b_2 & c_2 & ... & 0 \\ 0 & & \ddots & & 0 \\ 0 & ... & a_{n-1} & b_{n-1} & c_{n-1} \\ 0 & & ... & a_n & b_n\\ \end{bmatrix}.\] In some applications with periodic boundaries, a variant called the cyclic tridiagonal matrix occurs, \[A' = \begin{bmatrix} b_1 & c_1 & 0 & ... & \beta\\ a_2 & b_2 & c_2 & ... & 0 \\ 0 & & \ddots & & 0 \\ 0 & ... & a_{n-1} & b_{n-1} & c_{n-1} \\ \alpha & & ... & a_n & b_n\\ \end{bmatrix}.\] The efficient Thomas algorithm, and any optimized program we may have written to implement it, cannot be applied.

However, thanks to the Sherman-Morrison formula, we are still in business!

Following the lead from Numerical Recipes, we can set \[u = \begin{bmatrix} \gamma \\ 0 \\ \vdots \\ 0 \\ \alpha \end{bmatrix}, \quad \quad v = \begin{bmatrix}1 \\ 0 \\ \vdots \\ 0 \\ \beta/\gamma \end{bmatrix},\] and modify the elements of the tridiagonal part \(A\) of the cyclic tridiagonal matrix \(A'\), \(b_1^\prime = b_1 - \gamma\), and \(b_n^{\prime} = b_n - \alpha \beta/\gamma\).

Here \(\gamma\) is arbitrary, but is set to \(-b_1\) by default.

Thus \(A' = A + uv^T\).

To solve \(A' x = b\), in addition to matrix multiplication, we need to make two calls to the tridiagonal system. Thus, generally speaking, the asymptotic cost remains \(\mathcal{O}(n)\), but is about twice as costly as a standard tridiagonal system.