Is the Thomas algorithm the fastest way to solve a symmetric diagonally dominant sparse tridiagonal linear system

I am wondering if the Thomas algorithm is the fastest way (provably?) to solve a symmetric diagonally dominate sparse tridiagonal system in terms of algorithmic complexity (not looking for implementation packages like LAPACK etc). I know that both the Thomas algorithm and multigrid are $O(n)$ complexity, but maybe the constant factor for multigrid is less? It doesn't seem to me like multigrid could be faster but I am not positive.

Note: I am considering the case where the matrices are very large. Either direct or iterative methods are acceptable.

I believe comparing an iterative method (multigrid) to a direct/exact method (Thomas) in terms of exact operation count isn't really meaningful. IIRC, Thomas operation count is $8N$ for any tridiagonal system. The only time I can imagine multigrid conceivably beating that is for a trivial case of having a linear solution, and even then the cost of evaluating the residual at each level would be comparable to the cost of Thomas.

The $O(N)$ usefulness of multigrid lies in the fact that it's general for sparse matrices, and not restricted to tridiagonal systems.

• Thanks. I realize that iterative methods are not exact. I should have specified a very small tolerance (say 10^-15) and just treated that as being "exact" for comparison purposes. – James Apr 28 '14 at 16:32
• @user2697246 well, you asked about "provably" fastest. The exact convergence rate for multigrid (or any iterative scheme) is always going to depend on the solution itself and the starting guess - a linear solution will be effectively solved exactly in one step, whereas something more oscillatory will take more operations. Thomas has an exact, fixed operation count for all cases. Practically speaking, you're never going to beat Thomas for (serially) solving a tridiagonal system for a non-trivial case. – Aurelius Apr 28 '14 at 16:48
• @Aurelius Can the Thomas algorithm be parallelized? If not, that is one major advantage of multigrid! – Nick Alger Apr 28 '14 at 17:21
• @NickAlger No, the Thomas algorithm is strictly serial, and yes parallelization is a big advantage to multigrid (although for the specific case of a tridiagonal system I suspect communication latency would kill you.) There is a technique specific to tridiagonal systems called parallel cyclic reduction (PCR) which is $O(N log N)$, parallelizable by $N$, which I have used successfully on GPUs. – Aurelius Apr 28 '14 at 17:25
• One correction, the Thomas algorithm requires 8N operations, not 9N. Also, what do you mean by "multigrid ... having a linear solution"? All systems under consideration here are linear. – Doug Lipinski Apr 30 '14 at 13:23

The short answer is that the Thomas algorithm will be faster than any iterative scheme for almost all cases. The exception would perhaps be applying a single iteration of a very simple iterative scheme such as Gauss-Seidel, but this is highly unlikely to give an acceptable solution. Also, this is ignoring parallel processing concerns.

Multigrid is an especially poor choice in the case of a tri-diagonal matrix because although multigrid is $\mathcal O(n)$, the constant is quite large. In fact, multigrid does not even have an advantage over Gauss-Seidel until the matrices become quite large. This is due to the need for projection, prolongation, and relaxation operations for each multigrid level, each of which requires $\mathcal O(n)$ operations where n is the number of unknowns at that multigrid level.

Finally, this question is best addressed via operation counting. For the Thomas algorithm, a total of $5N$ multiplications and $3N$ additions are required for the solution. Iterative schemes require at least as many operations as matrix-vector multiplication and given a tri-diagonal matrix, each matrix-vector multiplication requires $3N-2$ multiplications and $2N-2$ additions. Therefore, even two applications of any (even the very simplest) iterative scheme will be more expensive than the Thomas algorithm.

• "Multigrid is an especially poor choice in the case of a tri-diagonal matrix because although multigrid is O(n), the constant is quite large." I think this too, but googling brought up a line in Trottenburg's Multigrid book claiming a constant of 0.1-0.2, stated without proof. I don't think I believe that. – Aurelius Apr 28 '14 at 15:38
• @Aurelius Interesting. That's clearly impossible in the general case since there are 3N entries in a tridiagonal matrix. If the cost is ~0.1*N, that means you never even operate on most of the entries. – Doug Lipinski Apr 28 '14 at 15:52
• Yes we're on the same page; simply evaluating a 3-point stencil requires 3N operations. I was just skimming so maybe I totally misinterpreted the statement, but you can see it for yourself in the google books excerpt. – Aurelius Apr 28 '14 at 16:40
• The full quote (pg 21) is "Efficiency in the practical sense means that the proportionality constants in this O(N) statement are small or moderate. This is indeed the case for multigrid: if designed well, the h-independed convergence factors can be made very small (in the range of 0.1-0.2 or even less) and the operation count per unknown per iteration step is also small." The 0.1-0.2 is refering to the residual reduction for each cycle of multigrid. The constant on the O(N) would be on the order of 1.5-2.0x the matrix multiply per cycle (with a total of a dozen or two cycles). – Godric Seer Apr 28 '14 at 18:15
• Ah, thanks @GodricSeer, that makes more sense. – Aurelius Apr 28 '14 at 18:20

Multigrid loops even on single core are vectorizable by the optimizer. So while operation counts can help, we should not forget that even in the serial world, processors have vector parallelism, and hence time-to-solution might not be exactly what we predict from the cost analysis.