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

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


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


10

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


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


8

In my experience, you always need (or better use) some form of preconditioning. The type and complexity of the precondition would vary depending on the task though. From Y. Saad, Iterative Methods for Sparse Linear Systems: Preconditioning is a key ingredient for the success of Krylov subspace methods... Lack of robustness is a widely recognized weakness of ...


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

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

Wolfgang Bangerth's answer already says almost everything, but another subtle detail is that GMRES/MINRES minimize the norm of the residual, i.e., $\|Ax_k-b\|$, while CG minimizes the (A-)norm of the error, i.e., $\|x_*-x_k\|_A$. In most cases what you really want to minimize is the error, and the residual serves as an imperfect proxy: recall that by the ...


5

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


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


5

In some sense it doesn't matter: In finite dimensions, all norms are equivalent, so if an algorithm conveniently has the property that proving convergence in one norm is easy, then that's what people will choose. In practice, one oftentimes wants to reduce the norm of the residual by a certain factor, and from a practical perspective, a certain reduction in ...


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


4

How GMRES method finds smallest singular value and the corresponding singular vectors of a matrix? It doesn't. GMRES solves linear systems. Your citation probably refers to other Krylov methods: restarted Arnoldi, Golub-Kahan bidiagonalization.


3

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


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


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

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

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


2

In case it may be of help to others, here's the result of the testing on a 30x30 influence matrix.It was made using non-restarted GMRES. GMRES without preconditioner (red curve) converges in 13 iterations, with a sudden convergence at iteration 12. GMRES with symmetric Gauss-Seidel (dashed curve) converges more smoothly, but in a little more iterations. ...


2

It seems like you have $N$ systems of linear equations that you need to solve. Let's say that you have $P$ cores available to you and unlimited memory. Now, you have two options: You can try solving $P$ systems in parallel, while each system is solved using 1 core. You can solve one system using $P$ cores, solving all $N$ systems one-by-one. The choice ...


2

Well, I can only propose a potential cause for why GMRES method fails for the problem you showed. I don't have enough reputation so I can't comment. Since GMRES is using Anorldi Process to generate a set of orthogonal vectors in the Krylov subspace, if the matrix itself is very large, that means the Q and H matrix produced could be very large. So if you ...


2

Taking $x_0 = 0$, we have that $x_1 \in <M^{-1}b>$. For the next iteration, we get $x_2 \in <M^{-1}b, M^{-1}AM^{-1}b>$. For the next iteration, we get $x_3 \in <M^{-1}b, (M^{-1}A)^2M^{-1}b>$. Continuing this argument, you will eventually find that $x_k \in \mathcal{K}_k(M^{-1}A,M^{-1}b)$. This space can also be described as $Mx_k \in \...


1

First some comments on why such matrices are hard for solvers. Notice that if $U$ is upper diagonal like [[. 2 . . .] [. . 2 . .] [. . . 2 .] [. . . . 2] [. . . . .]] then $(I - U)^{-1} = I + U + U^2 + ... U^{n-1} \ (U^n = 0)$. For this example, $(I - U)^{-1} =$ [[ 1 2 4 8 16] [ . 1 2 4 8] [ . . 1 2 4] [ . . . 1 2] [ . . . . 1]] ...


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