14

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


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


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


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

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


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

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


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


9

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


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

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


8

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


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


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


5

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


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


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

You might be familiar with the following paper already: http://link.springer.com/chapter/10.1007%2F978-3-642-22061-6_10 Problems which are highly indefinite and oscillatory are very difficult to design robust iterative methods for. The paper gives some suggestions which might be helpful to you though, many of them have been extended to the time-harmonic ...


4

In short, orthogonalization of the Krylov vectors occurs with respect to the operator, but not with respect to the preconditioner. Alright, so say we want to solve $Ax=b$ with preconditioner $B$. the preconditioned-CG iteration is basically: \begin{align*} \hat{v}_1=\tilde{v}_1 =& Bb\\ v_1 =& \tilde{v}_1 / c_1\\ \\ \hat{v}_i =& BAv_{i-1}\\ \...


4

If your matrices are large, why not use a library like ARPACK? The shift-and-invert mode of ARPACK will help you calculate the eigenvalues close to $\sigma$. There are interfaces to ARPACK for most high-level programming languages used for numerical computations (Fortran, C, Python, MATLAB, etc.). For example, a quick tutorial for Python can be found here ...


4

It can be even harmful. In Liesen/Strakos Krylov subspace methods principles and analysis (Chapter 5.8.3) it is reported that a nonzero initial x0 makes a GMRes iteration first remove unwanted components of x0 before it starts the approximation of the wanted solution. If a nonzero initial guess is used, it should be rescaled by the Hegedus trick; see the ...


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

Convergence of iterative methods are affected by the condition number of the matrix, which tends to increase as the mesh is refined. Benzi has done work showing that improved convergence can be obtained using ILU or approximate inverse preconditioners if one first performs permutations to maximize the diagonal entries of the matrix. Search for "...


3

Both ILU and diagonal scaling are not efficient preconditioners for "real" problems, i.e., if you let the size of your discrete problem become large. In other words, if you insist on timing GPUs vs CPUs with these preconditioners, you are only considering academic problems that are not really of interest -- if the problem is small enough that the use of ...


3

Since expokit seems to use a Krylov subspace method, usually (at least, the hope is that) the upper Hessenberg matrices are of small dimension, say $m \sim 100$. For matrices of these sizes, there should not be any major difference in computational time by using any method for dense matrix exponential computation. For example, 'expm' in MATLAB seems to use ...


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