Cost analysis of Thomas algorithm

Question

I am doing somework using the Thomas algoritm when is comes to solving tridiagonal matrix. However, I cant find any information regarding the computational cost of the Thomas algorithm (TA). I do know that the computational cost of TA is O(n), but I just dont know how to "prove" it or show it. My task is to compare TA to gaussian elmination with pivoting.

Is there anyway I can show that the computaional cost of TA is O(n)?

Answer 1

You can just count operations... Assuming you have a tridiagonal matrix with diagonal $b_1, \cdots, b_n$, sub-diagonal $a_2, \cdots, a_n$ and super-diagonal $c_1,\cdots, c_{n-1}$, the Thomas algorithm consists in computing $c', d'$ according to $$ c_1'=c_1/b_1, \quad c_i'=\frac{c_i}{b_i-a_ic_{i-1}'}, \,\, i=2,\cdots, n-1 $$

$$ d_1'=d_1/b_1, \quad d_i'=\frac{d_i-a_id_{i-1}'}{b_i-a_i c_{i+1}'}, \,\, i=2,\cdots, n $$

and then computing $x$ by back-substitution

$$ x_n = d_n',\quad x_i = d_i'-c_i'x_{i+1}, \,\, i = n-1, \cdots, 1 $$

So, counting operations you have:

• Sums/subtractions: $(n-2) + 2(n-1)+(n-1) = 4n-5$
• Multiplications/divisions: $(1+2(n-2)) + (1 + 3(n-1)) + (n-1) = 6n-6$

So, the total operation count is $10n -11$, clearly $O(n)$, as mentioned in the comments.

Answer 2

I find it difficult to count flops by scanning the mathematical statement of many algorithms. There is one method that significantly reduces the risk of mistakes.

1. Implement the algorithm and verify that the software returns accurate results.
2. Augment the software to also count the number of flops.
3. Analyze the software manually and formalize a conjecture regarding the flop count.
4. Verify experimentally that the conjecture agrees with the exact flop count.

We must complete Step 1 before we do anything else. Step 2 is straight forward if we examine the source code one line at a time. Once we have implemented the algorithm and examined the arithmetic operations is much simpler to complete Step 3. If we discern that the flop count is a polynomial of degree $k$, then we need no more than $k+1$ experiments to verify that the we have isolated the correct polynomial.

In the case of the Thomas algorithm for solving a tridiagonal linear system of dimension $m$ the exact flop count is $8m-7$.

Below follows a MATLAB implementation of the algorithm that adheres to the principles listed above. Residuals are explicitly calculated in order to catch grievous programming errors and compare against the theoretical expectations. The flop count refers to the flop count for the factorization and solve phases.

The factorization function has a flop count of $3m-3$. We observe that this count is consistent with the case of $m=1$ where the factorization is trivial and does not require any flops.

function [L, U, R, rres, flops, diff]=thomasFactor(A)

% Isolate dimension of the problem
m=size(A,1);

% Initialize flop count 
flops=0;

% Create a copy
B=A;

% Main loop has length m-1
for j=1:m-1
    % Compute multiplier
    B(j+1,j)=B(j+1,j)/B(j,j);
    % Update the trailing submatrix;
    B(j+1,j+1)=B(j+1,j+1)-B(j+1,j)*B(j,j+1);
    % Update flops
    flops=flops+3;
end

% Extract L and U
L=eye(m,m)+tril(B,-1); U=triu(B);

% Compute residual
R=A-L*U;

% Compute relative residual
rres=norm(R,'inf')/norm(A,'inf');
    
% Exact count based on analyzing the main loop.
count=3*m-3;

% Difference
diff=count-flops;

The solve function has a flop count of $5m-4$. We observe that this count is consistent with the case of $m=1$ where the solve requires a single division.

function [x, r, rres, flops, diff]=thomasSolve(A,L,U,b)

% Isolate the dimension of the problem
m=size(L,1);

% Initialize flop count
flops=0;

% Initialize solution
x=b;

% Forward sweep of length m-1
for j=2:m
    % Do the linear update
    x(j)=x(j)-L(j,j-1)*x(j-1);
    % Update flops
    flops=flops+2;
end

% Do the central division
x(m)=x(m)/U(m,m);

% Update flops 
flops=flops+1;

% Backward sweep of length m-1
for j=m-1:-1:1
    % Do the linear update
    x(j)=x(j)-U(j,j+1)*x(j+1);
    % Do the central division
    x(j)=x(j)/U(j,j);
    % Update flops
    flops=flops+3;
end

% Compute residual
r=b-A*x;

% Compute relative residual
rres=norm(r,'inf')/norm(b,'inf');

% Exact count based on analyzing the loops
count=5*m-4;

% Difference 
diff=count-flops;