# What's wrong with the **PCG and MINRES** in matlab?

Last week, I have learned the details of the robust iterative methods of PCG, MINRES, GMRES, which will converges to the exact solution $$x^*$$ of nonsingular system within $$N$$ steps for $$A\in \mathbb{R}^{N\times N}$$: $$Ax=b.$$ And the requirements for the matrix is: SPD matrix for PCG, symmetric matrix for MINRES, nonsymmetric matrix for GMRES. The reason is that these 3 above methods are optimal, i.e., their residual vectors are mutually orthogonal, so their residual vectors, e.g., $$r_0,r_1,...,r_N$$ must satisfy $$r_N=0$$, which produces the exact solution in the $$N$$-th step. So, I have done some examples in matlab (matlab 2018b, 8GB win10 OS environment) as follows:

clc;clear;
rng(0);%    to fix the rand number
n=10;
A=rand(n);  %   generate a nonsingular matrix
b=rand(n,1);
%   for general matrix
gmres(A,b);
%%  for SPD matrix
pcg(A*A',b);
minres(A*A',b);
gmres(A*A',b);


And the results are:

gmres converged at iteration 10 to a solution with relative residual 0.

pcg stopped at iteration 10 without converging to the desired tolerance 1e-06
because the maximum number of iterations was reached.
The iterate returned (number 10) has relative residual 0.02.

minres stopped at iteration 10 without converging to the desired tolerance 1e-06
because the maximum number of iterations was reached.
The iterate returned (number 10) has relative residual 0.016.

gmres converged at iteration 10 to a solution with relative residual 0.


It is so strange that gmres is normal for both nonsymmetric and SPD matrix but pcg and minres are abnormal. Because matrix $$AA^T$$ is SPD, then in theory, pcg and minres must converge in N steps, but they fail. Does it mean that we should use GMRES as much as possible instead of other methods regardless of the SPD matrix or else, because gmres command is so robust for nonsymmetric and SPD matrix? Is anything wrong? any suggestions?

• pcg converges in 11 iterations for me in octave – Nick Alger Dec 9 '19 at 3:04
• $r_N=0$ in exact arithmetic, but due to floating point precision and the resulting degradation of the conjugation of search directions, for example, this will rarely hold in practice. – cdipaolo Dec 9 '19 at 3:04
• @cdipaolo Seems strange that loss of orthogonality would be such a big issue for such a small matrix. When I try it, the answer at 10 iterations is horrible, but the answer at 11 iterations is accurate to machine precision. Maybe there is an off-by-1 thing here – Nick Alger Dec 9 '19 at 3:11
• @cdipaolo, thanks for your reply, I get it due to the loss of orthogonality of the residual vector, I have done the examples that pcg converges in 11 steps, and minres also 11 steps, but gmres still in 10 steps, so gmres is more useful for SPD matrix or general matrix, right? – sunshine Dec 9 '19 at 4:46
• I think it is really loss of orthogonality what's going on here. If you wish, you can check for yourself: code up the Arnoldi iteration for that matrix, and compute $Q^T Q - I$. You will see that its entries grow in magnitude as one gets away from the diagonal. – Federico Poloni Dec 9 '19 at 23:13

While very similar, each method is slightly different and you should definitely take this into account.

The GMRES method is the simplest, it will construct an orthogonal basis for $$\mathcal{K}_k(A,b)$$ and select an approximate solution that minimises the 2-norm of the residual.

The MINRES method is a variation of GMRES based on the fact that for a symmetric matrix, the vector that will become the new basis vector will already be orthogonal to all but the last 2 previous basis vectors. This leads to a large reduction in computation time since we don't need to do all that orthogonalisation. However, in floating point arithmetic, we might lose some of that orthogonality (this is also true for GMRES by the way, mgs is technically not stable).

The PCG method is similar but instead of being the symmetric variant of GMRES, it is a variant of FOM.

Now on to your problem. As already mentioned in the comments, all methods do converge after 11 steps. However, I've had a look at the residual vector and the internal code of matlab and an off-by-1 does not seem to be the problem. One of the causes might be that matrices constructed as $$AA^*$$ have a condition number that the square of the original condition number of $$A$$. This can lead to slower convergence and a worsening of the loss of orthogonality.

It might be something else too, I'm not entirely sure, but I hope my answer has made it clear that each method has its advantages and disadvantages. Concluding that GMRES is better in every way is a step too far.

• Thanks Prof. Steel for your clear and detailed reply, I will go on studying Krylov subspace methods. and thanks for your help. – sunshine Dec 11 '19 at 10:26
• You're welcome, but I'm not a professor ;) – Thijs Steel Dec 11 '19 at 12:42
• Though not a Prof., however, your understanding of Krylov subspace methods is very excellent. In the future, I will ask for many help about theories of Krylov subspace. Thanks again for your generosity. ^_^ – sunshine Dec 12 '19 at 0:41