Hot answers tagged

15

Interesting that this question came yesterday, since I just finished an implementation yesterday that does this. My Background Just to start of, let me know that while my education background is from scientific computing, all work I have done since graduating, including my current Ph.D. work, has been in computational electromagnetics. So, I guess our ...


11

There is one major difference between GMRES over FOM. It is also the reason why I would recommend GMRES over FOM. In exact arithmetic, the residuals obtained by GMRES form a decreasing sequence. You are certain that the GMRES residuals will not increase in the absence of rounding errors. Once the computed residual deviates from this simple pattern there is ...


9

GMRES is indeed one of the hardest to understand Krylov methods. As you correctly state, the algorithm computes in each step $m$ a new approximation $x^m$ to the solution of $Ax=b$ as a minimizer of $\|b-Ax\|_2$ over the Krylov space $K_m(A,b) = \mathrm{span}\{b,Ab,A^2b,\dots,A^{m-1}b\}$. (I'm assuming for simplicity that $x^0 =0$). Clearly, if $A\in\mathbb{...


9

Such a nested Krylov subspace method may work quite well in practice. It may be of interest for non symmetric linear systems for which restarted GMRES stagnates and unrestarted GMRES is too expensive or uses too much memory. Some literature: GMRESR: A family of nested GMRES methods, van der Vorst, Vuik Flexible inner-outer Krylov subspace methods, Simoncini,...


9

The MUMPS sparse direct solver can handle symmetric indefinite systems and is freely available (http://graal.ens-lyon.fr/MUMPS/). Ian Duff was one of the authors of both MUMPS and MA57 so the algorithms have many similarities. MUMPS was designed for distributed-memory parallel computers but it also works well on single-processor machines. If you link it ...


8

One thing that CG has in its favor is that it's not minimizing the discrete $l^2$ norm for its residual polynomials (what GMRES does). It's minimizing a matrix-induced norm instead, and very often this matrix-induced norm ends up being very close to the energy norm for discretizations of physical problems, and frequently this is a much more reasonable norm ...


7

The general problem that direct solvers are suffering from is the fill-in phenomenon, meaning that the inverse of a sparse matrix may be dense. This leads to huge memory requirements if the structure of the matrix is not "suitable". There are attempts to work around these issues, and MATLAB's default lu-function employs a few of them. See http://...


7

Performing $k$ steps of GMRES uses $O(n k^2)$ time and $O(n k)$ memory. In other words, the algorithm gets more and more expensive with each additional iteration. In theory, the algorithm terminates in $k \le n$ steps (ignoring round-off), thereby yielding a worst-case cubic $O(n^3)$ time and quadratic $O(n^2)$ memory figure. In practice, we manually ...


7

Unfortunately, convergence of GMRES does not have a clear dependence on the distribution of eigenvalues. It was proved by Greenbaum, Ptak and Strakos in 1996 that you can construct examples with an arbitrary spectrum and an arbitrary convergence history: that is, give me any $n$ nonzero complex numbers, and any decreasing sequence $\|r_k\|$, and I can ...


7

We can't use the criterion you show last in practice because it requires us to know what $\kappa(A)$ is. But computing the condition number is, in general, more expensive than solving a linear system. As a consequence, the criterion you show is not practical. That only leaves us with variants of the criterion $$ \frac{\|r_k\|}{\|r_0\|}=\frac{\|b-Ax_k\|}{\|b-...


6

The reason is that GMRES can only be used for solving linear equations, i.e. equations of the form $Ax=b$, where $A$ is some matrix and $x,b$ are vectors. What GMRES does, essentially, is it approximates multiplication by the matrix $A^{-1}$ using a matrix polynomial of $A$. In this case (I assume) $f(y^{n+1},t)$ is not necessarily linear in the vector $y^{...


6

This is all wasted effort if you are not using a preconditioner because no choice of parameters will lead to a robust or scalable algorithm, and the least-bad choice will be effected by preconditioning. After you have found a good preconditioner, you can look at the effect of restarts and decide on how to manage restarts in GMRES, as well as comparing to ...


5

Both PETSc and Trilinos have high quality implementations. Both of these libraries also have bindings to other languages.


5

It seems that you are referring to the fluxes that occur at the edge of your domain. Since these are typically known quantities enforced by your boundary conditions, there is no need to solve for those values. There are two common ways to deal with this: Include an equation for each (known) boundary quantity in the form $f_i = F(x_i)$ for $x_i\in\partial\...


5

Technically, the template Matlab code provided at Netlib already calculates the obtained residual at every iteration. It's just not recorded in the way that this is an output of the subroutine. After each iteration iter the residual r is calculated, so one just needs to record it. I will do it right after it is calculated in the residual vector resvec(iter)....


5

Estimates for GMRES convergence based on eigenvalue distribution often implicitly assume that the matrix is normal. Sometimes the convergence rate is still provable in an asymptotic sense in the non-normal case, but if the matrix is severely non-normal then the "pre-asymptotic" behavior will make such convergence rates never reachable in practice. Your ...


5

First of all, MATLAB's gmres assumes that the preconditioner you use is linear. This is important! Actually it is the main difference between FGMRES and GMRES. Right preconditioned GMRES and FGMRES are exactly the same if you use a linear preconditioner, however, FGMRES allows the use of non-linear preconditioners. What do I mean by a non-linear ...


4

It is indeed a tolerance on the relative residual $$ ||b-Ax_k|| / ||b||$$ In the webpage that you mention, you can find this information (very implicitly) [x,flag,relres] = gmres(A,b,...) also returns the relative residual norm(b-A*x)/norm(b). If flag is 0, relres <= tol. The third output, relres, is the relative residual of the preconditioned ...


4

I suspect there is in general not much difference between GMRES and CG for an SPD matrix. Let's say we are solving $ Ax = b $ with $ A $ symmetric positive definite and the starting guess $ x_0 = 0 $ and generating iterates with CG and GMRES, call them $ x_k^c $ and $ x_k^g $. Both iterative methods will be building $ x_k $ from the same Krylov space $ K_k ...


4

In addition to @jed-brown's answer, I can highly recommend the book Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics by Elman, Silvester and Wathen. The book covers finite elements, but the general strategies for Krylov subspace methods also hold for finite volumes, of course.


4

The intuition behind all Krylov's subspace methods is the following. Given a square matrix $A$ and a compatible vector $b$, there exists a unique monic polynomial $p$ such that $p(A)b = 0$. Mathematically, this is a consequence of the axiom of choice and Cayley's theorem which ensures that $q(A) = 0$ where $q$ is the characteristic polynomial of $A$. If $A$ ...


4

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 ...


3

You can often manipulate systems like this using a block LU factorization; the factors you mentioned then come into play. Preconditioners can then be designed by taking advantage of the block diagonal structure revealed by the factorization, with the blocks being $M_1$ and the Schur complement $M_2-S^TM_1^{-1}S$. An advantage of this approach is that you ...


3

One thing is that GMRES is never used wherever CG can be applied. I dont think it does make sense to compare these two. For SPD matrices, CG is definitely the winner because of the storage requirements and the reasons you mentioned above. A question that would be interesting is, to find an extension of CG, that is applicable to problems where CG can not be ...


2

Any iterative solver can beat direct methods only if the problem is sufficiently large (large, depends on several factors such as storage required, efficiency of implementation). And also any krylov method (for example GMRES) are good only if you use an appropriate preconditioning (in practice). If you are not using any preconditioning, krylov methods are of ...


2

Matlab '\' operator is highly efficient due to top notch programming. You may see some of the idea and how they used all possible short cuts in Timothy Davis's book. In matlab you can use gmres, which is well developed and stable. Probably minres, which theoretically should be ideal for your case, may not be that reliable (at least my experience says so). ...


2

From what I can see, on a normal machine, you would run out of memory if you are using a Dense matrix. The implementation you indicated in the question uses a dense matrix ( storing all the matrix entries ). Generally in most simulations, each element computes fluxes with its adjacent elements and the resulting matrix structure is sparse ( the number of non-...


2

I'm not sure I understand the question, but in GMRES you built a orthonormal base for the space $K_m \; \text{ for } m= 1 \dots N$. In this space the generic vector can be write as: $$ x = x_0 + V_m y$$ where: $x_0 \in \mathbb{R}^N $ $\text{dim}(V_m) = N \times m$ $y \in \mathbb{R}^m$, $\text{dim}(y) = m \times 1$ During the various steps the algorithm ...


2

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: Construct a sequence of approximations by solving a projected least squares problem on a nested sequence of subspaces $\mathcal{K}_m$. Choose the nested subspaces specifically as the ...


2

Please check this paper by Benzi et al. They address this issue and give corresponding references on p. 45. Shortcut: for the Stokes problem $A = \text{diag}(A_{11},A_{11},\dots,A_{11})$ is just a collection of discrete Laplace operators, so it is natural to approximate their inverses using multigrid. However, things get much more complicated for Oseen type ...


Only top voted, non community-wiki answers of a minimum length are eligible