The Thomas Algorithm

Suppose you need to solve $n$ equations in $n$ unknowns, but each equation involves only its immediate neighbors — the unknown just before and just after it. Write that as a matrix and you get a tridiagonal matrix: every entry is zero except on the main diagonal and the two adjacent diagonals.

Systems like this appear constantly in science and engineering. Discretize a heat equation, a beam deflection, or a fluid velocity profile on a grid and the resulting linear system is almost always tridiagonal. Solve thousands of such systems per second and you care deeply about cost.

Naive Gaussian elimination costs $O(n^{3})$ operations. But for a tridiagonal matrix almost all of those operations touch a zero and accomplish nothing. Lev Thomas (1949) noticed that the band structure lets you perform elimination in a single forward pass and recover the solution in a single backward pass — $O(n)$ total, with a tiny constant.

The Thomas algorithm is exact (no approximations), numerically stable for diagonally dominant systems, and simple enough to fit on one page. It is one of the most practically important linear-time algorithms in numerical computing.

The demo below shows a 5-equation tridiagonal system. The three diagonals are labeled a (lower), b (main), c (upper), and d is the right-hand side.

<!-- {{c_html_intro}} -->
<p class="hint">{{hint_para}}</p>
<div id="matrix-display" class="matrix-display"></div>
<div class="status" id="status">{{status_ready}}</div>
<div class="btns">
  <button id="btn-sweep" type="button">{{btn_sweep}}</button>
  <button id="btn-back" type="button" disabled>{{btn_back}}</button>
  <button id="btn-reset" type="button" class="ghost">{{btn_reset}}</button>
</div>

/* {{c_css_intro}} */
* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; font-size: 14px; }
.hint { font-size: .88rem; color: #444; margin: 0 0 .7rem; line-height: 1.5; }
.matrix-display { display: grid; grid-template-columns: auto auto auto auto auto; gap: 3px 6px; margin: .5rem 0 .6rem; width: fit-content; }
.col-label { font: 700 11px ui-monospace, monospace; color: #888; text-align: center; padding-bottom: 2px; }
.row-label { font: 700 11px ui-monospace, monospace; color: #888; display: flex; align-items: center; padding-right: 2px; }
.cell { width: 54px; height: 36px; display: flex; align-items: center; justify-content: center;
        font: 600 13px ui-monospace, monospace; border-radius: 6px; border: 1px solid #cdd9e3;
        background: #e8eef3; color: #1d3557; transition: background .25s, color .25s; }
.cell.zero { color: #bbb; background: #f4f6f8; }
.cell.active { background: #1d3557; color: #fff; border-color: #1d3557; }
.cell.changed { background: #d4edda; color: #0a7d33; border-color: #a8d5b5; }
.cell.solution { background: #fff3cd; color: #856404; border-color: #ffc107; }
.cell.rhs { background: #f0e8f8; color: #5a1d9e; border-color: #c4a5e8; }
.cell.rhs.changed { background: #e8d5f8; }
.cell.rhs.solution { background: #fff3cd; color: #856404; border-color: #ffc107; }
.status { font-size: .95rem; font-weight: 600; margin: .4rem 0; min-height: 1.5em; }
.status.ok { color: #0a7d33; }
.status.info { color: #1d3557; }
.status.done { color: #856404; }
.btns { display: flex; gap: .5rem; flex-wrap: wrap; margin-top: .3rem; }
button { font: 600 13px system-ui, sans-serif; padding: .4rem .85rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 8px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
button:disabled { opacity: .4; cursor: default; }

// Code not found

Click Forward sweep to watch the algorithm eliminate the lower diagonal one row at a time, updating the main diagonal and right-hand side in place. Click Back-substitute once the sweep is done to recover the solution $x$ from bottom to top. Reset starts over with the original system.

Notice that no row ever touches anything outside its immediate neighbors — that is the structural miracle that makes the whole thing $O(n)$.

How fast is the Thomas algorithm, really? Let us count precisely.

Forward sweep. For each row $i$ from $2$ to $n$, compute a multiplier $w = a_{i} / b_{i-1}$, then set $b_{i} \leftarrow b_{i} - w \cdot c_{i-1}$ and $d_{i} \leftarrow d_{i} - w \cdot d_{i-1}$. That is 4 arithmetic operations per row: $n - 1$ rows, so $4(n-1)$ operations.

Back substitution. Compute $x_{n} = d_{n} / b_{n}$, then for $i = n-1$ down to $1$ set $x_{i} = (d_{i} - c_{i} \cdot x_{i+1}) / b_{i}$. That is 3 operations per step: $3(n-1) + 1$ total.

Grand total: $\approx 7n$ floating-point operations — firmly $O(n)$, versus the $O(n^{3})$ of general Gaussian elimination and even the $O(n^{2})$ of banded LU with a wide band. The saving for $n = 10{,}000$ is a factor of roughly $n^{2}/7 \approx 14{,}000{,}000$.

Stability. The algorithm is numerically stable without pivoting when the matrix is strictly diagonally dominant — $|b_{i}| > |a_{i}| + |c_{i}|$ for every row — a condition automatically satisfied by most finite-difference discretizations. For non-dominant systems cyclic reduction or pivoting variants are used instead.

Space. Only the three diagonals and the right-hand side are stored: $O(n)$ memory, with no fill-in. Compare general LU, which can require $O(n^{2})$ storage for the factors.

The Thomas algorithm is a special case of dynamic programming: each step depends only on the immediately preceding step, so no information needs to be carried further back.

Any problem whose unknowns interact only with their nearest neighbors on a line produces a tridiagonal system — and that covers an enormous slice of applied computing:

• Finite-difference PDEs: discretizing the 1-D heat equation, wave equation, or Poisson equation with a standard stencil yields a tridiagonal system at every time step. The Thomas algorithm is the inner loop.
• Cubic spline interpolation: finding the second derivatives that make a natural cubic spline pass smoothly through $n$ data points requires solving a tridiagonal system once. Every graphics, CAD, and animation pipeline uses this.
• Implicit ODE integrators: methods like Crank–Nicolson and backward Euler convert an ordinary differential equation into a sequence of tridiagonal solves. Implicit methods are stable for stiff problems — and fast precisely because of Thomas.
• Sweep methods in 2-D: alternating-direction implicit (ADI) schemes split a 2-D PDE into a sequence of 1-D tridiagonal problems, each solved by one Thomas pass. This extends linear-time solving into higher dimensions.
• Signal processing: certain IIR filters and recursive smoothers have tridiagonal structure; the Thomas algorithm processes them in one sweep.

Whenever you see a simulation or solver described as "implicit" and "efficient," there is almost certainly a Thomas algorithm — or its block-tridiagonal generalization — running underneath. It connects directly to ideas in dynamic programming and to the general study of algorithm design.

The Thomas algorithm is a lesson in the power of structure. A general $n \times n$ linear system looks intractable at scale — $O(n^{3})$ operations, $O(n^{2})$ memory. But impose the rule that each unknown talks only to its two neighbors and almost all of that work vanishes. Two tidy passes — one forward, one backward — suffice, and the total cost is just $7n$ arithmetic operations.

The insight generalises: whenever a problem has the shape of a chain — each piece depending only on its immediate predecessor and successor — a linear-time sweep exists. That is the core idea behind dynamic programming, behind Viterbi decoding, and behind the Thomas algorithm itself.

Next time a physical model or a curve-fitting problem leads you to a tridiagonal system, you are not staring at a hard problem. You are holding one of the most efficient exact solvers ever devised.