Solve Tridiagonal matrix in for loop

Welcome!

I was looking for the fastest possible way to solve the triagonal equation. In my program, calling the part that solves a triagonal equation occurs many times because the values of vector B change every iteration of the for loop (the number of loop iterations reaches up even to several thousand times).

Initially, due to the lack of knowledge about other possible solutions, I used the Thomas algorithm available on the web. However, the operating times were far from satisfactory, so I wrote my own interpretation which significantly accelerated the work of the code.

While searching for even better ways to solve the problem, I came across the Linsolve function and the backslash solution, highly praised in the general forum for its speed.

Initially, I was enthusiastic about modifying the code in order to implement a simple linsolve function, but to my surprise, the code runtime increased several times compared to the Thomas version. That suprised me, so I made a short program that returns me the times of operation of individual solutions of the triagonal equation and compared the results for the selected methods:

load("test_sparse.mat");
Test = diag(D) + diag(E,1) + diag(C,-1);
max_iter = 1000;
%%
B_init = zeros(1,100);
[eta,h_init,r_init] = deal(B_init);
N = length(D);
[eta C_eta] = licz_eta(C,D,E,N);
tic;
for i =1:max_iter
wynik_macierzy = thomas(E,B,N,eta,h_init,r_init,C_eta)';
end
Time_thomas = toc;

%%

tic;
for i =1:max_iter
wynik_mldivide = mldivide(Test,B');
end
Time_mldivide = toc;

tic;
for i =1:max_iter
wynik_linsolve = linsolve(Test,B');
end
Time_linsolve = toc;
end

%%

Test_sparse = sparse(Test);
B_sparse = sparse(B);
tic;
for i =1:max_iter
wynik_backsl = Test/B';
%wynik_sparse(1) = wynik_sparse(1) + 25;
end
Time_backsl = toc;

For a single inclusion of the program, I got the following results:

Here it is very clear that the Thomas method is slightly slower than other solutions, so using them should bring benefits.

However, the situation changes when I set up more iterations of the loop. For 10,000 iterations:

Here Thomas wins with a big advantage, leaving the competition far behind.

Here comes my question. Why are methods that for a single iteration gave better time results with more loop iterations suddenly become much worse?

In addition, I noticed that all the functions discussed give the same final result, except for the backslash method, which looks simillar compared to the others but differs in the first line of the vector. I don't know where this difference comes from.

In the attachment all the files necessary for the code to work, if someone wanted to check it out.

Thank you very much in advance for all your time and help.

Matlab ver: R2020a

CPU: i5 6600K

[Christine Tobler]

Backslash will use a special-case solver for tridiagonal if the input matrix is sparse. It's not checking for this in a dense matrix, since it's more efficient to use a sparse matrix to store this case.

Use Test = sparse(Test) (more efficiently, use spdiags to create Test directly), and you should see improvement in the mldivide call. Sparse matrices are not supported in linsolve, but I don't think you need it: linsolve shouldn't be expected to be faster than mldivide unless additional information about the matrix structure is added.

[Mariusz Kijas]

I used sparse version of Test just like you recommend me. For now i calculate only runtime of function mldivide - rest of times such as creating sparse matrix is not counted here, of course in main program i will use spdiag.

Test_sparse = sparse(Test);
tic;
for i =1:max_iter
    wynik_mldivide_sparse = mldivide(Test_sparse,B');
end
Time_mldivide_sparse = toc;

The given change shows the following result (10,000 iterations) :

Indeed, replacing the matrix Test with its sparse version significantly improved the timework of the mldivide function. Unfortunately Thomas still produces better results.

[Christine Tobler]

Those timings make sense to me. Since mldivide works for general sparse matrices, it needs to check if the input matrix is tridiagonal on each step. Those

Also, the algorithm used in mldivide verifies if row-permutations are needed, to make sure the solution is sufficiently accurate. I don't think the Thomas algorithm does this - and if you know that your Test matrix is always positive definite, it's indeed not necessary to test this.