# Pcg method is imprecise

2 views (last 30 days)
Fedor on 7 Mar 2023
Commented: Fedor on 11 Mar 2023
I must be wrong in my code, so feel free to note where I've mistaken. I'm testing pcg method's precision through following scheme: generate random positive definitive matrix, and let it be linear system with the desired solution as (1, 1, 1, 1...) - to simplify testing. We can solve that linear system with division or simply state that solution should be (1, 1, 1, 1...). Then we solve system again with pcg method and test that all xi are nearby 1. And it doesn't work - pcg's solution does not meet needed tolerance(precision). Maybe I understand the meaning of tolerance incorrectly?
Code:
tol = 1e-8;
iterations = 150;
matrix_1_size = 50;
matrix_1 = sprand(matrix_1_size,matrix_1_size,.5);
matrix_1 = matrix_1'*matrix_1;
B_1 = sum(matrix_1, 2);
expected_solution_1 = matrix_1 \ B_1;
pcg_solution_1 = pcg(matrix_1, B_1, tol, iterations);
diff = abs(pcg_solution_1 - expected_solution_1);
tolerance = zeros(matrix_1_size, 1) + tol;
assert(all(diff <= tolerance));

Chunru on 8 Mar 2023
The pcg method use tol as a relative tolerance such that norm(residual) < tol*norm(rhs). So you need a different assert.
tol = 1e-8;
iterations = 150;
matrix_1_size = 50;
matrix_1 = sprand(matrix_1_size,matrix_1_size,.5);
matrix_1 = matrix_1'*matrix_1;
B_1 = sum(matrix_1, 2);
expected_solution_1 = matrix_1 \ B_1;
pcg_solution_1 = pcg(matrix_1, B_1, tol, iterations);
pcg converged at iteration 74 to a solution with relative residual 9.8e-09.
diff = abs(pcg_solution_1 - expected_solution_1);
norm(diff)
ans = 1.3296e-06
norm(pcg_solution_1-expected_solution_1)
ans = 1.3296e-06
norm(B_1)
ans = 772.9681
norm(pcg_solution_1-expected_solution_1)/norm(B_1) < tol
ans = logical
1
assert(norm(pcg_solution_1-expected_solution_1)/norm(B_1) < tol)
##### 1 CommentShow NoneHide None
Fedor on 11 Mar 2023
Thanks, that is the solution!

### Categories

Find more on Matrix Indexing in Help Center and File Exchange

R2022b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!