Why MATLAB chooses the Householder in its built-in function gmres.m?

Question

Recently, I have studied how to construct an orthonormal basis for Krylov subspace to solve $Ax=b$, where $A\in \mathbb{R}^{n\times n}$ is nonsingular. As we know, there are usually 4 ways to construct an $m$ dimensional Krylov subspace $\{v,Av,...,A^{m-1}v\}$ as follow:

1. standard Gram-Schmidt (G-S);
2. Modified Gram-Schmidt (MGS);
3. Householder reflection (House);
4. MGS with reorthogonalization (MGSR).

$$\begin{array}{|c|c|c|} \hline \textbf{Method} & \textbf{Work} & \textbf{Storage} \\ \hline \text{G-S} & 2m^2n & (m+1)n \\ \hline \text{MGS} & 2m^2n & (m+1)n \\ \hline \text{House} & 4m^2n-\frac{4}{3}m^3& (m+1)n-\frac{1}{2}m^2 \\ \hline \text{MGSR} & 4m^2n & (m+1)n \\ \hline \end{array}$$

My question is from the computational work, the MGSR is most about twice than G-S and MGS. From the numerical stability point, Householder is the most reliable method, but the computational work is very much. For us, if given any matrix $A$, how should we choose the best algorithm to write the method, e.g., gmres.m? Is there a criterion? By the way, I find my Matlab 2018b chooses the Householder method. Does this mean that in practice, Householder is proved the best efficient and stable method?

Answer 1

If you're using GMRES, typically you have a large stiff system. The extra work done for the householder algorithm is negligible compared to the expense of GMRES and the preconditionder. As such, we want the more numerically stable method to make sure that the system is more likely to converge. Especially since you choose GMRES because you want something that is guaranteed to converge, introduced numerical instability is suboptimal.