In this first and half week, as written in the first post, I studied the algorithms and also their codes in Octave. I recall that the algorithms are gmres, pcg, bicg, bicgstab, qmr, cgs and pcr.
The codes that I used that are not in Octave can be found at
https://bitbucket.org/CrisDorigo/socis16-octave-iter_meth
in the folder "First_week".
Octave codes
Studying the implementations I noticed problems in some of them:
- pcg - Trying to solve the simple problem 1*x = 1, with initial guess 1, there is the following output:
octave4.1-cli:29> x = pcg(1,1,[ ],[ ],[ ],[ ],1)
warning: division by zero
warning: called from
pcg at line 367 column 10
pcg: converged in 0 iterations. pcg: the initial residual norm was reduced NaN times.
x = 1
We can see that in 0 iteration there is a warning for a division by zero and that the residual is reduced NaN times. Then this output text needs to be fixed. - bicg - I tried to solve the small linear system A*x = b, with
A = [7 7 8 4 4; 7 0 6 3 4; 3 2 3 7 6; 6 0 9 7 7; 1 0 0 1 7];
b = sum(A,2);
(I obtained this matrix with floor(10*rand(5)) so it is a simple integer matrix, and b computed as the sum of the columns of A in such a way the exact solution is the vector x = [1; 1; 1; 1; 1])
With this configuration I discovered that the algorithm stops before the convergence, because of stagnation. Indeed the output is:
octave4.1-cli:39> [x, flag, relres, iter, resvec] = bicg(A,b)
bicg stopped at iteration 4 without converging to the desired tolerance 1.000000e-06
because the method stagnated.
The iterate returned (number 3) has relative residual 3.772464e-02
x =
1.67650
1.04764
0.47523
0.92726
0.97121
flag = 3
relres = 0.038052
iter = 4
resvec =
0.070581 0.042765 0.037725
As we can see the last residual (relres = 0.038052) is greater than the third (that is the third component of resvec, i.e. 0.037725). Then the algorithm detect a stagnation and stops.
I tried to solve the same linear system with Matlab, and I noticed that also its bicg has the fourth iteration with the residual greater than the third one, but it doesn't stop and at the fifth iteration it converges at the solution.
So I tried to comment the part of the stagnation (i.e. lines from 187 to 195 in the Octave's bicg) and tried to solve this linear system, obtaining:
octave4.1-cli:41> [x,flag,relres,iter,resvec] = bicg_1(A,b)
bicg stopped at iteration 5 without converging to the desired tolerance 1.000000e-06
because the maximum number of iterations was reached. The iterate returned (number 5) has relative residual 3.371209e-14
x =
1.0000
1.0000
1.0000
1.0000
1.0000
flag = 1
relres = 3.3712e-14
iter = 5
resvec =
0.070581 0.042765 0.037725 0.038052
Except that an error text that says that the maximum number of iterations was reached, we can see that at the fifth iteration the algorithm converges with tolerance smaller than the default one 1e-6.
The bicg with the "stagnation" commented is the file bigc_1.m, instead the matrix A and the test with bigc and bicg_1, can be found in example_bicg.m, in the folder First_week . - cgs - I studied the algorithm from the Saad's book, and then I tried to study the codes, but I noticed that its Octave implementation was not taken from this book. Since in the help there is not the reference book I can't compare this implementation with the reference.
I also tried to solve the same linear system discussed in the bicg section, and I noticed that the algorithm does not converge. The obtained output is:
octave4.1-cli:56> [x,flag,relres,iter,resvec] = cgs(A,b)
x =
0.9799298
-0.0060874
0.7904128
0.6945865
1.8885042
flag = 1
relres = 0.17700
iter = 5
resvec =
1.000000
0.070581
0.129415
0.158692
0.159693
0.177003
We can see that the solution given by the algorithm is far from the exact solution, also the flag equal to 1 indicates that the algorithm reached the maximum number of iterations without convergence.
Then I tried to make a simple and basic implementation from scratch using the formulation in the Saad book, I tried to solve the linear system with this implementation and this implementation converges. The output is:
octave4.1-cli:58> [x,resvec] = cgs_1(A,b)
x =
1.0000
1.0000
1.0000
1.0000
1.0000
resvec =
5.1604e+01
3.1724e+00
3.3616e+00
1.4778e+00
3.0564e+00
3.6586e-10
To be sure that is not a pathological case I tried to solve with the Octave implementation and mine the linear system A*x = b, with
A = [6 6 8 7 2; 5 4 1 0 3; 4 1 2 6 3; 6 9 3 8 3; 3 5 6 9 4];
b = sum(A,2);
obtaining similar results.
Then, after a discussion with the mentors, we decided that maybe it is easier to implement this method from scratch, instead to search the problems, due to the fact that there is no reference for its implementation.
My simple and basic cgs implementation is the file cgs_1.m, and the test with the second matrix can be found in the file example_cgs.m, in the folder First_week. - bicgstab - I noticed that in Matlab are displayed also the "half iteration" with the "half residual", but not in Octave. Then talking with the mentors we decided to display them also in Octave.
Matrix-vector product in the algorithms
Since all of the algorithms of the project are solvers for linear systems, the matrix-vector product is essential in all of their implementation, but since the matrix A and the (eventual) preconditioners can be passed as matrices or as function handle, it is not easy to decide how to compute all these matrix-vector products.
For example in the gmres algorithm the strategy is to transform A, M1, M2 into function handle, i.e. (tanking A as example) to create the function handle
Ax = @(x) A*x; if A is a matrix
Ax = A; if A is a function handle
and then every time that there is a matrix-vector product, for example A*v, it is simply computed as Ax(v) or feval(Ax, v).
Instead, in the pcg implementation, every time that there is a matrix-vector product it is checked if A is a matrix or a function handle. For example if we need to compute A*v then there are the following lines of codes:
if (isnumeric(A))
w = A*v;
else
w = feval(A,v);
endif
(the other algorithms use one of these two strategies, or similar variations).
Then it is natural the question: which is the best strategy?
To answer, I simulate these strategy and the simple matrix-vector product, then I computed the times, using sparse matrices of different dimensions (n = 100,
n = 200, n = 500, n= 1000, n = 2000) in this way
A = sprand(n, n, rand(1));
Ax = @(x) A*x;
w = rand(length(A),1);
# Estimation of the simply A*w
tic
for j =1:100
u = A*w;
endfor
t1 = toc;
# Estimation of feval(Ax, w)
tic
for j = 1:100
u = feval (Ax, w);
endfor
t2 = toc;
# Estimation of the if block with A as matrix
tic
for j = 1:100
if (isnumeric (A))
u = A*w;
else
u = feval(A,w);
endif
endfor
t3 = toc;
# Estimation of the if block with Ax the function handle
tic
for j = 1:100
if (isnumeric (Ax))
u = Ax*w;
else
u = feval(Ax,w);
endif
endfor
t4 = toc;
obtaining (for example, since I used random matrices) as results:
n = 100
t1 = 0.0038838
t2 = 0.0068600
t3 = 0.0049579
t4 = 0.0076430
n = 200
t1 = 0.0051990
t2 = 0.0083339
t3 = 0.0062981
t4 = 0.0091600
n = 500
t1 = 0.0042169
t2 = 0.0071800
t3 = 0.0051830
t4 = 0.0081360
n = 1000
t1 = 0.055708
t2 = 0.059194
t3 = 0.056590
t4 = 0.060054
n = 2000
t1 = 0.15554
t2 = 0.16691
t3 = 0.16327
t4 = 0.16677
This test with different matrices can be found in the file quicktest.m in the folder First_week.
We can see that if we consider n<=1000, the best strategy is to use the simple matix-vector product (indeed t1 is the less value). Instead if we consider the two strategies of gmres and pcg, we notice that the "if block" used in pcg takes more time than the simple matrix-vector product and the evaluation of the function handle.
Instead the case n = 2000 is not a good example to see the differences of these approaches, probably because the matrix is too big and requires a lot of operations.
Then, after a discussion with the mentors, we decided that the "best strategy" is to write in every algorithm two nested function "ad hoc" to the type of the data, and if A, M1 and M2 are not of the same type, to transform all of them in function handle. To be clear, for example, as sketch of how it will be:
function [x, flag, relres, iter, resvec] = pcg (A, b, tol, maxit, M1, M2, x0)
if A, M1, M2 are matrices
[x, flag, relres, iter, resvec] = pcg_matrix (A, b, tol, maxit, M1, M2, x0)
if A, M1, M2 are function handle
[x, flag, relres, iter, resvec] = pcg_function (A, b, tol, maxit, M1, M2, x0)
if A, M1, M2 are not of the same type
[x, flag, relres, iter, resvec] = pcg_mixed (A, b, tol, maxit, M1, M2, x0)
endfunction
Priority of the algorithms
The last part of this post is about the chosen priority of the improvements of the algorithms for the next weeks.
The mentors suggested me to improve the algorithms in the following order:
- pcg
- gmres
- bicgstab
- pcr
- qmr
- cgs
- bicg
Then follows pcr because it is the only method already implemented in Octave that solve linear system when the input matrix A is Hermitian (bicgstab and gmres work for general matrices, instead pcg works for SPD matrices).
And in the last follows qmr, cgs and bicg because they are the less used methods, with cgs above bicg because it has to be written from scratch.
If there are something that is not clear or if there are some suggestions/advices, comment below or contact me via email at cristiano.dorigo@hotmail.it
