2
$\begingroup$

The GMRES method is for solving the linear system $Ax=b$. Given an initial guess $x_0$ and the corresponding residual $r_0:=b-Ax_0$, we have the Krylov subspace $$\mathcal{K}_m:=\mathop{span}\{r_0,Ar_0,A^2r_0,\ldots,A^{m-1}r_0\}.$$ The GMRES method tries to find an approximate solution $x_m$ with the correction $d_m:=x_m-x_0$ from $\mathcal{K}_m$ such that $\|r_0-Ad_m\|_2$ is minimized. The GMRES method uses the Arnoldi process to get a decomposition $$AV_m=V_{m+1}H_m$$ with $V_{m+1}$ an orthonormal basis of $\mathcal{K}_{m+1}$ and $H_m$ an $(m+1)$-by-$m$ Hessenberg matrix. Let $d_m=V_my$ with $y\in\mathbb{R}^m$. Then, we have $$\|r_0-Ad_m\|_2=\|V_{m+1}^Tr_0-H_my\|_2$$ which leads to a least square problem that is smaller and easier to solve.

My question is that why not use directly the original basis $X_m:=[r_0,Ar_0,A^2r_0,\ldots,A^{m-1}r_0]$. Let $d_m=X_mz$ with $z\in\mathbb{R}^m$. Then, we do QR decomposition directly for $AX_m$ so that $$AX_m=QR$$ with $R$ an upper triangular matrix. Then, using this decomposition we can still minimize $\|r_0-Ad_m\|_2$ easily.

$\endgroup$
  • 1
    $\begingroup$ The first step in any practical QR decomposition is to transform the matrix to upper Hessenberg form, which the Arnoldi process gives you for free. Also, the practical implementation of the GMRES algorithm makes use of the QR decomposition of $H_m$ (actually, the extended matrix) to compute the decomposition of $H_{m+1}$ with a single Givens rotation and postpones actually solving the least squares system until the residual would be small enough. All this exploits the properties of the Arnoldi process. $\endgroup$ – Christian Clason Jul 3 '17 at 6:50
  • $\begingroup$ "The first step in any practical QR decomposition is to transform the matrix to upper Hessenberg form". But Householder can directly do QR without the first step. "Also, the practical implementation of the GMRES algorithm ...". With my proposal, I think I can also postpone solving the least square problem-- just doing the Householder for the new column each time. $\endgroup$ – Hui Zhang Jul 3 '17 at 7:35
  • $\begingroup$ Another thing which I thought was obvious but is worth pointing out is that $H_{m}$ is an $(m+1)\times m$ matrix, while $AX_m$ is an $n\times m$ matrix. In practice, $m\ll n$ (think $m=10$ versus $n=10^4$). This also has bearing on your comment: Of course Householder can do it, but Givens is much more efficient in this case. $\endgroup$ – Christian Clason Jul 3 '17 at 8:11
  • $\begingroup$ But note that $H_m$ was obtained by the Arnoldi process which costs about the same as the Householder for $AX_m$. $\endgroup$ – Hui Zhang Jul 3 '17 at 8:38
  • 3
    $\begingroup$ Well, if you copy all the neat tricks of GMRES (column-wise update of Hessenberg form and QR factorization, delayed solution of triangular system), you end up with an algorithm with at best cosmetic differences, so what you're proposing is in the end not very different. The standard approach is likely more stable, though, since you orthogonalize each new Krylov basis vector as soon as you compute it instead of first forming $A^{m-1}r_0$ (which converges to the dominant eigenvector of $A$ for $m\to\infty$). Fundamentally, it's the difference between Gram-Schmidt and modified Gram-Schmidt. $\endgroup$ – Christian Clason Jul 3 '17 at 11:30
2
$\begingroup$

Since nobody else has answered this question (which is a good question, by the way!) so far, I'm consolidating my comments here. Briefly, GMRES boils down to three key ideas:

  1. Construct a sequence of approximations by solving a projected least squares problem on a nested sequence of subspaces $\mathcal{K}_m$.

  2. Choose the nested subspaces specifically as the Krylov spaces $\mathcal{K}_m=\mathrm{span} \{r^0,Ar^0,\dots,A^{m-1}r^0\}$.

  3. Solve the least squares problem efficiently using transformation to Hessenberg form and QR decomposition, exploiting the fact that the subspaces are nested and you only need to update a single column of the decomposition in each iteration.

Assuming you keep the general structure (in particular, the transformation to Hessenberg form and the column-wise update of $AX_m$ and its QR factorization), the only difference between GMRES and your approach is the following:

  • In your approach, you
    1. build the matrix $$AX_m = A[r^0, Ar^0,\dots,A^{m-1}r^0] = [Ar^0, A^2r^0,\dots,A^{m}r^0],$$
    2. transform it to Hessenberg form via Householder transformations, and
    3. compute (or update) its QR decomposition with Givens rotations.

Note that the QR decomposition implicitly constructs (as the columns of $Q$) an orthonormal basis of the space spanned by the columns of $AX_m$.

  • In comparison, GMRES
    1. builds (through the Arnoldi process) an orthonormal basis of $X_m$ and a transformation (via the same basis!) of $AX_m$ to upper Hessenberg form, and
    2. computes (or updates) its QR decomposition with Givens rotations.

So the only difference is that GMRES avoids the formation of the matrix $AX_m$ and instead directly proceeds to its Hessenberg form. In effect, you are immediately orthogonalizing every new Krylov vector as soon as it is formed rather than waiting until the end. This is relevant because as in the power iteration, the vectors $A^mr^0$ converge to the dominant eigenvector of $A$, so that the matrix $AX_m$ becomes increasingly ill-conditioned. If you wish, you can say that GMRES is the modified Gram-Schmidt (which is more stable) to your standard Gram-Schmidt.

(And if you now say, "well, but I can do the same: orthogonalize $A^mr^0$ before adding it as a column to $AX_m$": congratulations, you've just reinvented GMRES.)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.