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


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

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


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

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


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

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


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

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

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


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

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

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.


4

As far as I know, there are no such methods in LAPACK. Since LAPACK is the linear algebra package, no nonlinear solvers are included. However, you can use the underlying BLAS for implementing iterative methods. For nonlinear solvers using Jacobians (e.g., Newton's method), the matrix factorizations of LAPACK may come in handy. You may want to have a look at ...


3

This question is pretty well discussed in literature. However, there are lots of questions concerning multigrid on SciComp, so I decided to compose more or less detailed answer. I. When multigrid DOES work fine as a stand–alone solver MG works fine for rather simple elliptic problems. In fact, both MG and PCG (with MG cycles as inner iterations) usually ...


3

As mentioned in Geoff Oxberry's more complete answer, it should be noted that PETSc includes TSSUNDIALS, an interface to SUNDIALS. If you configure PETSc with the --download-sundials option (see python2 ./configure --help | grep -A 2 sundials for other related options, such as using an existing SUNDIALS library), then you can use at least some of the ...


3

In general, you can do more with PETSc. SUNDIALS is a collection of ODE solvers (in CVODE, Adams-Bashforth and BDF methods; in ARKODE, ARKIMEX methods) and DAE solvers (IDAS implements a BDF method) with sensitivity capabilities (the CVODES and IDAS variants), and a nonlinear solver (KINSOL). There are a few Krylov solvers in there (at least GMRES and CG) ...


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


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

In order to understand these results, you need to know how minimization and projection problems are connected. Namely, Let $\mathbb V$ be a subspace of $\mathbb C^n$ and take $y \in \mathbb C^n$; then $\hat y = \text{argmin}_{x \in \mathbb V} || y - x ||$ iff $(\hat y, x) = (y, x)$ for all $x \in \mathbb V$. Here $(.,.)$ is a scalar product (not ...


3

I doubt I can explain this better than the author, but I'll give it a shot. Let's say that $r_0 = \sum \alpha_i x_i$, with $x_i$ an eigenvector with eigenvalue $\lambda_i$. We can then write the vectors in the basis as $A^kr_0 = \sum \lambda_i^k\alpha_i x_i$. If all eigenvalues are distinct, $A^kr_0$ will converge to the eigenvector with the largest (in ...


3

Your view is correct, Matlab does store the zeros. As pointed out by @rchilton1980, this particular non-optimization that you are pointing out here is not too harmful, since the bulk of the storage in Krylov methods is the matrix V, not H. But that is just an instance of a general phenomenon. That is just Matlab's choice of a tradeoff between simplicity and ...


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