I need to solve a tridiagonal system (positive definite, diagonally dominant) $Ax = b$ in a time stepping loop. $A \in \mathbb{R}^{N \times N}$ remains constant but $b$ changes during each time step.

I was thinking which of the following would be faster :

  1. Compute and store $\mathbf{z_i}$ that solves $A\mathbf{z_i} = \mathbf{e_i} \; \forall 1\leq i\leq N$. Then any $b \in \mathbb{R}^{N}$ is a linear combination of $\mathbf{z_i}$ and $x = b_i \mathbf{z_1} + \dotsb b_N \mathbf{z_N}$

  2. Naively solve $A x = b$ for each given $b$ by tridiagonal matrix algorithm.

I thought (1) will save the time taken for computing the tridiagonal solves since the $z$ essentially gives the inverse of $A$. However, I calculated that one solve for $\mathbf{e_i}$ takes $8N$ flops, so that is $8N \times N$ to compute all the $\mathbf{z_i}$. Then for each given $b$, it takes $N^2$ flops to compute $x$. If I take $n$ time-steps, that is a total of $8N^2 + nN^2$.

In contrast (2) needs $n\times 8N$.

Am I correct ? Where is (2) gaining the advantage ?

  • 2
    $\begingroup$ Since $A$ remains constant, the usual preference is to store $A$ in its factored form $LL^T$, since the back-solve and forward-solve are $O(N)$ complexity. $\endgroup$ – hardmath Apr 2 '16 at 15:42
  • $\begingroup$ That $N^2$ may be missing a coefficient 2 in front of it, but your analysis looks correct to me. As usual, computing inverses explicitly is the wrong choice -- not at all surprising: the same happens for dense matrices. $\endgroup$ – Federico Poloni Apr 2 '16 at 15:47

The most efficient way would likely be to perform one step of Gauss-Elimination on A to eliminate the sub-diagonal and then store the inverse of the new main diagonal, the new modified super-diagonal and the constants to preform forward substitution on b.

This way each rhs will require one forward substitution to modify b followed by one back-solve. The forward operations consist of N additions and N multiplications and the backward operations consist of 2N multiplications and N additions. This adds up to 5N operations and no divisions (which are slower).

| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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