26

Initial advice Always run with -ksp_converged_reason -ksp_monitor_true_residual when trying to learn why a method is not converging. Make the problem size and number of processes as small as possible to demonstrate the failure. You often gain insight by determining what small problems exhibit the behavior that is causing your method to break down and the ...


21

In general, all Krylov methods essentially seek a polynomial that is small when evaluated on the spectrum of the matrix. In particular, the $n$th residual of a Krylov method (with zero initial guess) can be written in the form $$ r_n = P_n (A) b $$ where $P_n$ is some monic polynomial of degree $n$ . If $A$ is diagonalizable, with $A=V\Lambda V^{-1}$, we ...


19

Introduce the vector $y:=-A^{-1}Gx$ and solve the large coupled system $Ay+Gx=0$, $G^Ty=-b$ for $(y,x)$ simultaneously, using an iterative method. If $A$ is symmetric (as seems likely though you don't state it explicitly) then the system is symmetric (but indefinite, though quasidefinite if $A$ is positive definite), which might help you to choose an ...


18

Yes, you can, but Krylov methods generally do not have great smoothing properties. This is because they target the whole spectrum in an adaptive way that minimizes the residual or a suitable norm of the error. This will generally include some low frequency (long wavelength) modes that the coarse grids would have handled fine. Krylov smoothers also make the ...


17

On norms As an addendum to Reid.Atcheson's answer, I would like to clarify some issues regarding norms. At the $n^{\mathrm{th}}$ iteration, GMRES finds the polynomial $P_n$ that minimizes the $2$-norm of the residual $$r_n = A x_n - b = \big(P_n(A) - 1 \big)b - b = P_n(A) b .$$ Suppose $A$ is SPD, so $A$ induces a norm and so does $A^{-1}$. Then $$\begin{...


16

My advice to students is to try a direct solver in these cases. The reason is that there are two classes of reasons why a solver may not converge: (i) the matrix is wrong, or (ii) there is a problem with the solver/preconditioner. Direct solvers almost always yield something that you can compare with the solution you expect, so if the answer of the direct ...


13

Warning Solving saddle point problems involves a lot more choices than definite problems, and there are a lot more things that can go wrong. Use monitors for all levels to debug convergence, to sure that null spaces are handled correctly when auxiliary operators are singular (usually just a constant null space), and to ensure that preconditioners are ...


13

You can use additive $$ P_a^{-1} x = (B^T B)^{-1} x + (C^T C)^{-1} x, $$ multiplicative $$ P_m^{-1} x = (B^T B)^{-1} x + (C^T C)^{-1} \Big(x - A (B^T B)^{-1} x \Big), $$ or symmetric multiplicative. Methods of this class are available in PETSc using PCCOMPOSITE in PETSc. For example, petsc/src/ksp/ksp/examples/tutorials$ ./ex2 -m 100 -n 100 -...


12

To perform reasonably, polynomial preconditioners need fairly accurate spectral estimates. For ill-conditioned elliptic problems the smallest eigenvalues are usually separated such that methods like Chebyshev are far from optimal. The most interesting property of polynomial methods is that they do not require any inner products. It's actually quite popular ...


12

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

I am extremely surprised that there is no mention of conditioning or the shape of the spectrum in your discussion, as it will be the decisive property in whether or not iterative methods can beat dense methods. As an extreme example, suppose that your dense matrix is some small perturbation of the identity matrix. Then most iterative methods will converge ...


11

Iterative methods in a nutshell: Stationary methods are in essence fixed point iterations: To solve $Ax=b$, you pick an invertible matrix $C$ and find a fixed point of $$ x = x + Cb- CAx $$ This converges by Banach's fixed point theorem if $\|I-CA\|<1$. The various methods then correspond to a specific choice of $C$ (e.g., for Jacobi iteration, $C=D^{-1}...


11

Practical experience shows that trying to get good initial iterates has little value. For example, in the context of solving partial differential equations, if you take the solution from one mesh, interpolate it onto a finer mesh, and use that as the starting guess for something like a CG iteration to solve the same problem on the finer mesh, then it turns ...


10

This is called "structurally symmetric". It simplifies graph traversal, such as occurs when setting up aggregates in algebraic multigrid, but doesn't offer much structure to improve convergence rates. Note that all common PDE discretizations have this property so this is still a huge class of matrices including many instances for which no truly good ...


10

Computing the determinant of a sparse matrix is typically as expensive as a direct solve, and I am skeptical that CG would be of much help in computing it. It would be possible to run CG for $n$ iterations (where $A$ is $n \times n$) in order to generate information for the entire spectrum of $A$, and to then compute the determinant as the product of the ...


10

Some references on rounding error analysis of Krylov methods: Gutknecht, Martin H., and Zdenvek Strakos. "Accuracy of two three-term and three two-term recurrences for Krylov space solvers." SIAM Journal on Matrix Analysis and Applications 22.1 (2000): 213-229. Paige, Christopher C., and Zdenvek Strakos. "Residual and backward error bounds in minimum ...


10

In some cases, (F)MG provides an algorithm with optimal properties. For instance, properly tuned FMG can solve some elliptic problems in a small number of "work units", where a work unit is defined to be the computational effort required to express the problem itself - in this case the operations to form the residual $b-Ax$ on the finest grid. This is such ...


8

In most cases I know of where Krylov methods are used for dense problems, the operator is a low-rank perturbation of the identity (obtained by discretizing a continuum operator which is a compact perturbation of the identity). Such operators appear frequently in boundary integral equations as discretized Fredholm integral operators of the second kind. These ...


8

Your arguments apply naturally to the unpreconditioned case. The reason that I don't recommend pinning is that it confuses norms and preconditioning. If you know the size of a typical diagonal value, you can scale the trivial equation for the pinned node so that norms become reasonable again. To see the consequence on preconditioning, we have to distinguish ...


8

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


7

If your Krylov subspace is based on powers of $A$, convergence will be delayed by a number of iterations at most the rank of the correction. If it is based on powers of $A^TA$ then at most twice this number. I haven't seen such a statement in the literature. But to see the validity in the first case, it is sufficient to show that the $k$th Krylov space of ...


7

Following Arnold's reply, there is something you can do to simplify the problem. Specifically, rewrite the system as $Ay+Gx=0, G^Ty=-b$. Then note that from the statement that $G$ is tall and narrow and each row has only one 1 and zeros otherwise, then the statement $G^Ty=-b$ means that a subset of the elements of $y$ have a fixed value, namely the elements ...


7

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


6

For many partial differential equations arising in nature, particularly with strong nonlinearities or anisotropies, the choice of an appropriate preconditioner can have a large effect on whether the iterative method converges rapidly, slowly, or not at all. Examples of problems that are known to have fast and effective preconditioners include strongly ...


6

Others have already noted this but I think it's still worth pointing out that the determinant is not a useful quantity almost always when your matrix is large. The problem is that large matrices are most often approximations to things that are even larger dimensional (statistical samples of large populations, finite dimensional approximations to infinite ...


6

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


5

I think, this one Krylov Subspace Methods in Finite Precision: A Unified Approach, Jens-Peter M. Zemke is also worth reading.


5

CG was originally abandoned because of this loss of accuracy that you observed. It was revived when its use as an approximate solver was appreciated. So I don't think you want to revisit this problem. And preconditioning is very useful, diagonal preconditioning at the very least. Iterative methods are very problem dependent so I would recommend finding ...


5

The answer depends somewhat on the discretization; for example, some boundary integral discretizations result in very well-conditioned matrices, for which a Krylov solver works just fine without having to introduce the multi-level machinery of multigrid. For ill-conditioned matrices arising from finite difference or finite element methods, Krylov solvers ...


5

I must admit I never actually checked all the details myself, but I think that's a sketch of the general idea. The $k$th iterate $x_k$ produced by Richardson iteration lies in the Krylov subspace $K_k(A,b)$. The $k$th iterate $x_k$ produced by a Krylov method typically minimizes some objective function inside that same Krylov space, hence it is "better" ...


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