View source for Tridiagonal matrix algorithm

{{Short description|Improved reduction for specific matrices}}
In [[Numerical linear algebra|numerical linear algebra]], the '''tridiagonal matrix algorithm''', also known as the '''Thomas algorithm''' (named after Llewellyn Thomas), is a simplified form of [[Gaussian elimination]] that can be used to solve [[Tridiagonal matrix|tridiagonal systems of equations]]. A tridiagonal system for ''n'' unknowns may be written as

:<math>a_i x_{i-1}  + b_i x_i  + c_i x_{i+1}  = d_i,</math>
where <math>a_1 = 0</math> and <math>c_n = 0</math>.

:<math>
\begin{bmatrix}
   b_1 & c_1 &        &        &  0      \\
   a_2 & b_2 & c_2    &        &         \\
       & a_3 & b_3    & \ddots &         \\
       &     & \ddots & \ddots & c_{n-1} \\
   0   &     &        & a_n    & b_n
\end{bmatrix}
\begin{bmatrix}
   x_1    \\
   x_2    \\
   x_3    \\
   \vdots \\
   x_n
\end{bmatrix}
=
\begin{bmatrix}
   d_1    \\
   d_2    \\
   d_3    \\
   \vdots \\
   d_n
\end{bmatrix}
.
</math>

For such systems, the solution can be obtained in <math>O(n)</math> operations instead of <math>O(n^3)</math> required by [[Gaussian elimination]]. A first sweep eliminates the <math>a_i</math>'s, and then an (abbreviated) backward substitution produces the solution. Examples of such matrices commonly arise from the discretization of 1D Poisson equation and natural cubic [[Spline interpolation|spline interpolation]].

Thomas' algorithm is not [[Numerical stability|stable]] in general, but is so in several special cases, such as when the matrix is diagonally dominant (either by rows or columns) or symmetric positive definite;<ref name="Niyogi2006">{{cite book |author=Pradip Niyogi |title=Introduction to Computational Fluid Dynamics |year=2006 |publisher=Pearson Education India |isbn=978-81-7758-764-7 |page=76}}</ref><ref name="Datta2010"/> for a more precise characterization of stability of Thomas' algorithm, see Higham Theorem 9.12.<ref name="Higham2002">{{cite book |author=Nicholas J. Higham |title=Accuracy and Stability of Numerical Algorithms: Second Edition |year=2002 |publisher=SIAM |isbn=978-0-89871-802-7 |page=175}}</ref> If stability is required in the general case, [[Gaussian elimination]] with partial pivoting (GEPP) is recommended instead.<ref name="Datta2010">{{cite book |author=Biswa Nath Datta |title=Numerical Linear Algebra and Applications, Second Edition |year=2010 |publisher=SIAM |isbn=978-0-89871-765-5 |page=162}}</ref>

==Method==
The forward sweep consists of the computation of new coefficients as follows, denoting the new coefficients with primes:

: <math>c'_i =
 \begin{cases}
  \cfrac{c_i}{b_i},                  & i = 1, \\
  \cfrac{c_i}{b_i - a_i c'_{i - 1}}, & i = 2, 3, \dots, n - 1
 \end{cases}
</math>

and

: <math>d'_i =
 \begin{cases}
  \cfrac{d_i}{b_i},                                   & i = 1, \\
  \cfrac{d_i - a_i d'_{i - 1}}{b_i - a_i c'_{i - 1}}, & i = 2, 3, \dots, n.
 \end{cases}
</math>

The solution is then obtained by back substitution:
: <math>x_n = d'_n,</math>
: <math>x_i = d'_i - c'_i x_{i + 1}, \quad i = n - 1, n - 2, \ldots, 1.</math>

The method above does not modify the original coefficient vectors, but must also keep track of the new coefficients. If the coefficient vectors may be modified, then an algorithm with less bookkeeping is:

For <math>i = 2, 3, \dots, n,</math> do
: <math>w = \cfrac{a_i}{b_{i-1}},</math>
: <math>b_i := b_i-w c_{i-1},</math>
: <math>d_i := d_i-w d_{i-1},</math>

followed by the back substitution

: <math>x_n = \cfrac{d_n}{b_n},</math>
: <math>x_i = \cfrac{d_i - c_i x_{i+1}}{b_i} \quad \text{for } i = n - 1, n - 2, \dots, 1.</math>

The implementation as a [[C (programming language)|C]] function, which uses scratch space to avoid modifying its inputs for a-c, allowing them to be reused:
<source lang="c">
void thomas(const int X, double x[restrict X],
            const double a[restrict X], const double b[restrict X],
            const double c[restrict X], double scratch[restrict X]) {
    /*
     solves Ax = d, where A is a tridiagonal matrix consisting of vectors a, b, c
     X = number of equations
     x[] = initially contains the input v, and returns x. indexed from [0, ..., X - 1]
     a[] = subdiagonal, indexed from [1, ..., X - 1]
     b[] = main diagonal, indexed from [0, ..., X - 1]
     c[] = superdiagonal, indexed from [0, ..., X - 2]
     scratch[] = scratch space of length X, provided by caller, allowing a, b, c to be const
     not performed in this example: manual expensive common subexpression elimination
     */
    scratch[0] = c[0] / b[0];
    x[0] = x[0] / b[0];

    /* loop from 1 to X - 1 inclusive */
    for (int ix = 1; ix < X; ix++) {
        scratch[ix] = c[ix] / (b[ix] - a[ix] * scratch[ix - 1]);
        x[ix] = (x[ix] - a[ix] * x[ix - 1]) / (b[ix] - a[ix] * scratch[ix - 1]);
    }

    /* loop from X - 2 to 0 inclusive */
    for (int ix = X - 2; ix >= 0; ix--)
        x[ix] -= scratch[ix] * x[ix + 1];
}
</source>

==Derivation==
The derivation of the tridiagonal matrix algorithm is a special case of [[Gaussian elimination]].

Suppose that the unknowns are <math>x_1,\ldots, x_n</math>, and that the equations to be solved are:

:<math>\begin{alignat}{4}
&               && b_1 x_1 && + c_1 x_2       && = d_1 \\
& a_i x_{i - 1} && + b_i x_i && + c_i x_{i + 1} && = d_i \,, \quad i = 2, \ldots, n - 1 \\
& a_n x_{n - 1} && + b_n x_n &&                 && = d_n \,.
\end{alignat}
</math>

Consider modifying the second (<math>i = 2</math>) equation with the first equation as follows:

:<math>
(\mbox{equation 2}) \cdot b_1 - (\mbox{equation 1}) \cdot a_2
</math>

which would give:

:<math>
(b_2 b_1 - c_1 a_2) x_2  + c_2 b_1 x_3 = d_2 b_1 - d_1 a_2.
</math>

Note that <math>x_1</math> has been eliminated from the second equation. Using a similar tactic with the '''modified''' second equation on the third equation yields:

:<math>
(b_3 (b_2 b_1 - c_1 a_2) - c_2 b_1 a_3 )x_3 + c_3 (b_2 b_1 - c_1 a_2) x_4
 = d_3 (b_2 b_1 - c_1 a_2) - (d_2 b_1 - d_1 a_2) a_3.
\,</math>

This time <math>x_2</math> was eliminated. If this procedure is repeated until the <math>n^{th}</math> row; the (modified) <math>n^{th}</math> equation will involve only one unknown, <math>x_n</math>. This may be solved for and then used to solve the <math>(n - 1)^{th}</math> equation, and so on until all of the unknowns are solved for.

Clearly, the coefficients on the modified equations get more and more complicated if stated explicitly. By examining the procedure, the modified coefficients (notated with tildes) may instead be defined recursively:

:<math>\tilde a_i = 0\,</math>

:<math>\tilde b_1 = b_1\,</math>
:<math>\tilde b_i = b_i \tilde b_{i - 1} - \tilde c_{i - 1} a_i\,</math>

:<math>\tilde c_1 = c_1\,</math>
:<math>\tilde c_i = c_i \tilde b_{i - 1}\,</math>

:<math>\tilde d_1 = d_1\,</math>
:<math>\tilde d_i = d_i \tilde b_{i - 1} - \tilde d_{i - 1} a_i.\,</math>

To further hasten the solution process, <math>\tilde b_i</math> may be divided out (if there's no division by zero risk), the newer modified coefficients, each notated with a prime, will be:

:<math>a'_i = 0\,</math>

:<math>b'_i = 1\,</math>

:<math>c'_1 = \frac{c_1}{b_1}\,</math>
:<math>c'_i = \frac{c_i}{b_i - c'_{i - 1} a_i}\,</math>

:<math>d'_1 = \frac{d_1}{b_1}\,</math>
:<math>d'_i = \frac{d_i - d'_{i - 1} a_i}{b_i - c'_{i - 1} a_i}.\,</math>

This gives the following system with the same unknowns and coefficients defined in terms of the original ones above:

:<math>\begin{array}{lcl}
x_i + c'_i x_{i + 1} = d'_i \qquad &;& \ i = 1, \ldots, n - 1 \\
x_n = d'_n \qquad &;& \ i = n. \\
\end{array}
\,</math>

The last equation involves only one unknown. Solving it in turn reduces the next last equation to one unknown, so that this backward substitution can be used to find all of the unknowns:

:<math>x_n = d'_n\,</math>
:<math>x_i = d'_i - c'_i x_{i + 1} \qquad ; \ i = n - 1, n - 2, \ldots, 1.</math>

==Variants==
In some situations, particularly those involving [[Periodic boundary conditions|periodic boundary conditions]], a slightly perturbed form of the tridiagonal system may need to be solved:

:<math>
\begin{alignat}{4}
& a_1 x_{n}     && + b_1 x_1 && + c_1 x_2       && = d_1 \\
& a_i x_{i - 1} && + b_i x_i && + c_i x_{i + 1} && = d_i \,, \quad i = 2, \ldots, n - 1 \\
& a_n x_{n - 1} && + b_n x_n && + c_n x_1       && = d_n \,.
\end{alignat}
</math>

In this case, we can make use of the [[Sherman–Morrison formula]] to avoid the additional operations of Gaussian elimination and still use the Thomas algorithm. The method requires solving a modified non-cyclic version of the system for both the input and a sparse corrective vector, and then combining the solutions. This can be done efficiently if both solutions are computed at once, as the forward portion of the pure tridiagonal matrix algorithm can be shared.

If we indicate by:
<math display="block">
A=\begin{bmatrix}
   b_1 & c_1 &        &        &  a_1      \\
   a_2 & b_2 & c_2    &        &         \\
       & a_3 & b_3    & \ddots &         \\
       &     & \ddots & \ddots & c_{n-1} \\
   c_n   &     &        & a_n    & b_n
\end{bmatrix}
,x=
\begin{bmatrix}
   x_1    \\
   x_2    \\
   x_3    \\
   \vdots \\
   x_n
\end{bmatrix}
,d=
\begin{bmatrix}
   d_1    \\
   d_2    \\
   d_3    \\
   \vdots \\
   d_n
\end{bmatrix}

</math>

Then the system to be solved is: <math display="block">Ax = d</math>

In this case the coefficients <math> a_1 </math> and <math> c_n </math> are, generally speaking, non-zero, so their presence does not allow to apply the Thomas algorithm directly. We can therefore consider <math> B\in\mathbb{R}^{n\times n} </math> and <math> u,v\in\mathbb{R}^{n} </math> as following:
<math display="block">
B=\begin{bmatrix}
   b_1-\gamma & c_1 &        &        &  0      \\
   a_2 & b_2 & c_2    &        &         \\
       & a_3 & b_3    & \ddots &         \\
       &     & \ddots & \ddots & c_{n-1} \\
   0   &     &        & a_n    & b_n-\frac{c_na_1}{\gamma}
\end{bmatrix},
u=\begin{bmatrix}
   \gamma    \\
   0    \\
   0    \\
   \vdots \\
   c_n
\end{bmatrix}
,v=
\begin{bmatrix}
   1    \\
   0    \\
   0    \\
   \vdots \\
   a_1/\gamma
\end{bmatrix}
.
</math>
Where <math> \gamma\in\mathbb{R} </math> is a parameter to be chosen. The matrix <math> A </math> can be reconstructed as <math> A = B + u v^\mathsf{T}</math>. The solution is then obtained in the following way:<ref>{{Cite journal |last1=Batista |first1=Milan |last2=Ibrahim Karawia |first2=Abdel Rahman A. |date=2009 |title=The use of the Sherman–Morrison–Woodbury formula to solve cyclic block tri-diagonal and cyclic block penta-diagonal linear systems of equations |url=http://dx.doi.org/10.1016/j.amc.2009.01.003 |journal=Applied Mathematics and Computation |volume=210 |issue=2 |pages=558–563 |doi=10.1016/j.amc.2009.01.003 |issn=0096-3003}}</ref> first we solve two tridiagonal systems of equations applying the Thomas algorithm:
<math display="block">By=d \qquad  \qquad Bq=u </math> 

Then we reconstruct the solution <math>x </math> using the [[Sherman–Morrison formula|Shermann-Morrison formula]]:
<math display="block">\begin{align}
x &=A^{-1}d
=(B+uv^T)^{-1}d
=B^{-1}d-\frac{B^{-1}uv^TB^{-1}}{1+v^TB^{-1}u}d
=y-\frac{qv^Ty}{1+v^Tq}
\end{align} </math>


The implementation as a [[C (programming language)|C]] function, which uses scratch space to avoid modifying its inputs for a-c, allowing them to be reused:
<source lang="c">
void cyclic_thomas(const int X, double x[restrict X], const double a[restrict X], const double b[restrict X], const double c[restrict X], double cmod[restrict X], double u[restrict X]) {
    /*
     solves Ax = v, where A is a cyclic tridiagonal matrix consisting of vectors a, b, c
     X = number of equations
     x[] = initially contains the input v, and returns x. indexed from [0, ..., X - 1]
     a[] = subdiagonal, regularly indexed from [1, ..., X - 1], a[0] is lower left corner
     b[] = main diagonal, indexed from [0, ..., X - 1]
     c[] = superdiagonal, regularly indexed from [0, ..., X - 2], c[X - 1] is upper right
     cmod[], u[] = scratch vectors each of length X
     */

    /* lower left and upper right corners of the cyclic tridiagonal system respectively */
    const double alpha = a[0];
    const double beta = c[X - 1];

    /* arbitrary, but chosen such that division by zero is avoided */
    const double gamma = -b[0];

    cmod[0] = alpha / (b[0] - gamma);
    u[0] = gamma / (b[0] - gamma);
    x[0] /= (b[0] - gamma);

    /* loop from 1 to X - 2 inclusive */
    for (int ix = 1; ix + 1 < X; ix++) {
        const double m = 1.0 / (b[ix] - a[ix] * cmod[ix - 1]);
        cmod[ix] = c[ix] * m;
        u[ix] = (0.0f  - a[ix] * u[ix - 1]) * m;
        x[ix] = (x[ix] - a[ix] * x[ix - 1]) * m;
    }

    /* handle X - 1 */
    const double m = 1.0 / (b[X - 1] - alpha * beta / gamma - beta * cmod[X - 2]);
    u[X - 1] = (alpha    - a[X - 1] * u[X - 2]) * m;
    x[X - 1] = (x[X - 1] - a[X - 1] * x[X - 2]) * m;

    /* loop from X - 2 to 0 inclusive */
    for (int ix = X - 2; ix >= 0; ix--) {
        u[ix] -= cmod[ix] * u[ix + 1];
        x[ix] -= cmod[ix] * x[ix + 1];
    }

    const double fact = (x[0] + x[X - 1] * beta / gamma) / (1.0 + u[0] + u[X - 1] * beta / gamma);

    /* loop from 0 to X - 1 inclusive */
    for (int ix = 0; ix < X; ix++)
        x[ix] -= fact * u[ix];
}
</source>

There is also another way to solve the slightly perturbed form of the tridiagonal system considered above.<ref>{{Citation |last1=Ryaben’kii |first1=Victor S. |title=Introduction |date=2006-11-02 |url=http://dx.doi.org/10.1201/9781420011166-1 |work=A Theoretical Introduction to Numerical Analysis |pages=1–19 |publisher=Chapman and Hall/CRC |isbn=978-0-429-14339-7 |access-date=2022-05-25 |last2=Tsynkov |first2=Semyon V.|doi=10.1201/9781420011166-1 }}</ref> Let us consider two auxiliary linear systems of dimension <math>(n-1) \times (n-1) </math>:
<math display="block">
\begin{align}
 \qquad \ \ \ \ \ b_2 u_{2}  + c_2 u_3 &= d_2 \\
 a_3 u_2 + b_3 u_3 + c_3 u_4  &= d_3 \\ 
 a_i u_{i-1} + b_i u_i + c_i u_{i+1} &= d_i\\
 \dots \\
 a_n u_{n - 1}+ b_n u_n\qquad &= d_n \,.
\end{align}
\quad i = 4, \ldots, n - 1 \qquad \qquad


\begin{align}
 \qquad \ \ \ \ \ b_2 v_{2}  + c_2 v_3 &= -a_2 \\
 a_3 v_2 + b_3 v_3 + c_3 v_4  &= 0 \\ 
 a_i u_{i-1} + b_i u_i + c_i u_{i+1} &= 0\\
 \dots \\
 a_n v_{n - 1}+ b_n v_n\qquad &= -c_n \,.
\end{align}
\quad i = 4, \ldots, n - 1
</math>

For convenience, we additionally define <math> u_1 = 0</math> and <math>v_1 = 1</math>. We can now find the solutions <math>\{u_2,u_3\dots,u_n \} </math> and <math>\{v_2,v_3\dots,v_n \}</math> applying Thomas algorithm to the two auxiliary tridiagonal system.

The solution <math>\{ x_1,x_2\dots,x_n \}</math> can be then represented in the form:
<math display="block">x_i = u_i + x_1 v_i \qquad i=1,2,\dots, n</math>

Indeed, multiplying each equation of the second auxiliary system by <math>x_1</math>, adding with the corresponding equation of the first auxiliary system and using the representation <math>x_i = u_i + x_1 v_i</math>, we immediately see that equations number <math>2</math> through <math>n</math> of the original system are satisfied; it only remains to satisfy equation number <math>1</math>. To do so, consider formula for <math>i=2</math> and <math>i=n</math> and substitute <math>x_2 = u_2 + x_1 v_2</math>and <math>x_n = u_n + x_1 v_n</math> into the first equation of the original system. This yields one scalar equation for <math>x_1</math>:
<math display="block"> b_1x_1+c_1(u_2+x_1v_2)+a_1(u_n+x_1v_n) = d_1</math>

As such, we find:
<math display="block"> x_1 = \frac{d_1-a_1u_n-c_1u_2}{b_1+a_1v_n+c_1v_2}</math>

The implementation as a [[C (programming language)|C]] function, which uses scratch space to avoid modifying its inputs for a-c, allowing them to be reused:
<source lang="c">
void cyclic_thomas(const int X, double x[restrict X], const double a[restrict X], const double b[restrict X], const double c[restrict X], double cmod[restrict X], double v[restrict X]) {
    /* first solve a system of length X - 1 for two right hand sides, ignoring ix == 0 */
    cmod[1] = c[1] / b[1];
    v[1] = -a[1] / b[1];
    x[1] = x[1] / b[1];

    /* loop from 2 to X - 1 inclusive */
    for (int ix = 2; ix < X - 1; ix++) {
        const double m = 1.0 / (b[ix] - a[ix] * cmod[ix - 1]);
        cmod[ix] = c[ix] * m;
        v[ix] = (0.0f  - a[ix] * v[ix - 1]) * m;
        x[ix] = (x[ix] - a[ix] * x[ix - 1]) * m;
    }

    /* handle X - 1 */
    const double m = 1.0 / (b[X - 1] - a[X - 1] * cmod[X - 2]);
    cmod[X - 1] = c[X - 1] * m;
    v[X - 1] = (-c[0]    - a[X - 1] * v[X - 2]) * m;
    x[X - 1] = (x[X - 1] - a[X - 1] * x[X - 2]) * m;

    /* loop from X - 2 to 1 inclusive */
    for (int ix = X - 2; ix >= 1; ix--) {
        v[ix] -= cmod[ix] * v[ix + 1];
        x[ix] -= cmod[ix] * x[ix + 1];
    }

    x[0] = (x[0] - a[0] * x[X - 1] - c[0] * x[1]) / (b[0] + a[0] * v[X - 1] + c[0] * v[1]);

    /* loop from 1 to X - 1 inclusive */
    for (int ix = 1; ix < X; ix++)
        x[ix] += x[0] * v[ix];
}
</source>

In both cases the auxiliary systems to be solved are genuinely tri-diagonal, so the overall computational complexity of solving system <math> Ax = d </math> remains linear with the respect to the dimension of the system <math> n </math>, that is <math>O(n)</math> arithmetic operations.

In other situations, the system of equations may be '''block tridiagonal''' (see [[Block matrix|block matrix]]), with smaller submatrices arranged as the individual elements in the above matrix system (e.g., the 2D Poisson problem). Simplified forms of Gaussian elimination have been developed for these situations.<ref name="Quarteroni2007">{{cite book| last1=Quarteroni|first1=Alfio |last2=Sacco|first2=Riccardo | last3=Saleri|first3=Fausto | year=2007|title=Numerical Mathematics|publisher= Springer, New York|isbn= 978-3-540-34658-6 | chapter=Section 3.8}}</ref>

The textbook ''Numerical Mathematics'' by Alfio Quarteroni, Sacco and Saleri, lists a modified version of the algorithm which avoids some of the divisions (using instead multiplications), which is beneficial on some computer architectures.

Parallel tridiagonal solvers have been published for many vector and parallel architectures, including GPUs<ref name="ChangHwu14">{{cite conference |last1=Chang|first1=L.-W.|last2=Hwu|first2=W,-M.|title=A guide for implementing tridiagonal solvers on GPUs|book-title=Numerical Computations with GPUs | publisher=Springer|isbn=978-3-319-06548-9|editor=V. Kidratenko|year=2014}}</ref><ref name="Venetis">{{cite journal | last1=Venetis|first1=I.E. | last2=Kouris|first2=A.|last3=Sobczyk|first3=A.| last4=Gallopoulos|first4=E.|last5=Sameh|first5=A.|title=A direct tridiagonal solver based on Givens rotations for GPU architectures|journal=Parallel Computing|volume=49|pages=101–116|year=2015|doi=10.1016/j.parco.2015.03.008}}</ref>

For an extensive treatment of parallel tridiagonal and block tridiagonal solvers see <ref name="GPS.16">{{cite book| last1=Gallopoulos|first1=E. |last2=Philippe|first2=B.|last3=Sameh|first3=A.H. | year=2016|title=Parallelism in Matrix Computations|publisher= Springer|isbn= 978-94-017-7188-7|chapter=Chapter 5}}</ref>

==References==
<references />
* {{cite book |author1=Conte, S. D. |author2=de Boor, C. |year=1972 |title=Elementary Numerical Analysis |url=https://archive.org/details/elementarynumeri00samu |url-access=registration |publisher=McGraw-Hill, New York |isbn=0070124469}}
* {{Cite book |last1=Press |first1=W. H. |last2=Teukolsky |first2=S. A. |last3=Vetterling |first3=W. T. |last4=Flannery |first4=B. P. |year=2007 |title=Numerical Recipes: The Art of Scientific Computing |edition=3rd |publisher=Cambridge University Press |location=New York |isbn=978-0-521-88068-8 |chapter=Section 2.4 |chapter-url=http://apps.nrbook.com/empanel/index.html?pg=56}}

{{Numerical linear algebra}}

{{DEFAULTSORT:Tridiagonal Matrix Algorithm}}
[[Category:Numerical linear algebra]]

{{Sourceattribution|Tridiagonal matrix algorithm}}