14
$\begingroup$

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.

Improve this question
$\endgroup$
1
  • $\begingroup$ How can you apply a multigrid method to a tridiagonal system ? $\endgroup$
    – Yves Daoust
    Commented Oct 12, 2021 at 19:45

3 Answers 3

Reset to default
12
$\begingroup$

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.

Improve this answer
$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$
    – James
    Commented Apr 28, 2014 at 16:32
  • $\begingroup$ @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. $\endgroup$
    – Aurelius
    Commented Apr 28, 2014 at 16:48
  • $\begingroup$ @Aurelius Can the Thomas algorithm be parallelized? If not, that is one major advantage of multigrid! $\endgroup$
    – Nick Alger
    Commented Apr 28, 2014 at 17:21
  • 3
    $\begingroup$ @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. $\endgroup$
    – Aurelius
    Commented Apr 28, 2014 at 17:25
  • $\begingroup$ 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. $\endgroup$
    – Doug Lipinski
    Commented Apr 30, 2014 at 13:23
11
$\begingroup$

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.

Improve this answer
$\endgroup$
4
  • $\begingroup$ "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. $\endgroup$
    – Aurelius
    Commented Apr 28, 2014 at 15:38
  • 1
    $\begingroup$ @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. $\endgroup$
    – Doug Lipinski
    Commented Apr 28, 2014 at 15:52
  • $\begingroup$ 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. $\endgroup$
    – Aurelius
    Commented Apr 28, 2014 at 16:40
  • 5
    $\begingroup$ 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). $\endgroup$
    – Godric Seer
    Commented Apr 28, 2014 at 18:15
0
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.