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

Question

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?

Answer 1

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.