15

This question is very closely related to (and possibly a duplicate of) Is variable scaling essential when solving some PDE problems numerically?. There are still good practical reasons to nondimensionalize equations, if possible: It reduces the number of independent parameters for parametric studies (which was one of the original reasons for ...


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


15

Have a look at Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods (Barrett et al.). You can find it here. Here's why I'm recommending this over other references: the "flowchart of iterative methods" in appendix D (last page) covers both linear solvers and preconditioners, it is short (100 pages or so), does not go into too ...


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

Though it is a relatively rare situation when you actually have to calculate an inverse of the matrix, not all techniques were created equally. I would use the term badly-conditioned instead of ill-conditioned. For badly conditioned matrices, you might opt in the SVD-route to calculate the inverse: $$ A=U\Sigma V^H \implies A^{-1}=V\Sigma^{-1}U^H. $$ If ...


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


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

I'm going to answer my own question since the following method seems to be very effective. I'm making it an answer so people can upvote or downvote it independently of the question if they think it is good or bad. Answer: use randomized matrix probing applied to the diagonal of the matrix. Let $A$ be the operator we wish to find the diagonal of, and let $\...


9

The preconditioner for the FDM method that corresponds to the one you outline for the FEM (i.e., the Sylvester-Wathen approach) will still contain the Schur complement of the FDM matrix. The Schur complement of the FDM matrix will, in general, have the same kind of structure it has for the FEM, i.e., it will be spectrally equivalent to the identity operator ...


6

Depending on the problem structure, you can solve the ill-conditioned Augmented Lagrangian system directly. For example, BDDC/FETI-DP can solve almost-incompressible elasticity in primal form with a convergence rate independent of the Poisson ratio (piecewise constant on subdomains, but with arbitrary jumps). Similarly, multigrid methods that exactly ...


6

No, Jacobi only ever corrects relative scales. It does nothing for "smooth" ill-conditioning, such as the $\kappa(A) \in O(h^{-2})$ asymptotics for second order elliptic problems. If you are using a Krylov method, the global scale is automatically corrected, but with a stationary iteration, the (constant) scaling is needed somehow (could just be in the ...


6

This is something like a truncated SVD or eigenvector expansion of your solution. If you take $$x_m = \sum_{j=1}^m \frac{q_j\cdot b}{\lambda_j}q_j$$ with $m=n$, this is the exact solution to $Ax=b$. If you take instead the $m<n$ eigenvalues in the sum instead, you get an approximation. The problem is that it's a pretty poor approximation; since your ...


6

I don't know for sure whether there is an existing name for this method, but @jessechan's suggestion of "truncated eigenvector expansion" sounds perfectly fine to me (and most people would understand it). Your question had a second part embedded, namely why is this not what everyone does to get a good initial guess for iterative solvers? The answer to this ...


6

I don't particularly care for the notation $M^{-1}$ precisely because of the confusion you find yourself in. I (and others) simply call the preconditioner $P$. The point, however, is that for iterative solvers, you generally never need a matrix or its inverse explicitly. All you need to be able to do is apply it to a vector. So, let's say you want to solve $...


6

You should stick with GMRES, it is the only method that is essentially guaranteed to get a solution here. The real problem appears to be you need a better preconditioner. You could try sticking with LU but adding a diagonal mass matrix to the system with a constant multiplier that decreases with the linear residual of the system. This would allow you to get ...


5

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


5

Maybe not a preconditioning strategy in the traditional sense, but deflation could be useful in this case. In gmres(A) for instance, you can use the eigenpairs of the hessenberg projection H to form ritz vectors that are good estimates for eigenvectors of A. You use that to deflate your residual upon a restart, and give speedups over traditional restarted ...


5

I believe you're switching around the definitions of "inner" and "outer" stepping (this nomenclature gets worse in segregated schemes where you have at least 3 separate iterative schemes going on.) The core of what you're missing is that "pseudo-time" isn't a time at all; it's just a convenient way to achieve an iterative solution. Take a 1st order backward ...


5

You won't find any black-box scalable solutions because $S$ is typically dense and thus cannot be formed. If your problem comes from a mature research area, there might be experience in the literature demonstrating how to approximate Schur complements. One common technique is to use approximate commutator arguments. There are a quite a few papers on this ...


5

Based on your comments, it looks like your main difficulty is how having an LU factorisation for a given matrix $M$ makes it easy to find a vector $\mathbf{y}$ such that $M\mathbf{y}=\mathbf{b}$ for a given $\mathbf{b}$. Hopefully you're comfortable that if we have an easy method to do that, then this bit of the problem is solved, since we can treat $M^{-1}$ ...


5

Preconditioning and iterative solvers are cool, but did you try to solve your problem with some kind of sparse direct solver? If not, try it first. State–of–the–art preconditioning techniques for this kind of problems are multigrid techniques. In many cases multigrid allows you to solve your (elliptic) problem in $O(n)$ time. If you don’t want to go ...


5

If $\mathbf L$ and $\mathbf U$ give an approximate factorization of $\mathbf A$, you wouldn't want to use $\mathbf P = \mathbf L\cdot \mathbf U$ as a preconditioner (that's approximately $\mathbf A$), you'd want to use $\mathbf P = \mathbf U^{-1}\cdot \mathbf L^{-1}$ (that's approximately $\mathbf A^{-1}$). Even then, you would not want to multiply them ...


5

Starting with $AM^{-1} \mathbf{y} = \mathbf{b}$, where $\mathbf{y} \equiv M \mathbf{x}$, we can manipulate to $\mathbf{y} - (I-AM^{-1})\mathbf{y} = \mathbf{b}$ Replacing one instance of $\mathbf{y}$ by $\mathbf{y}^k$ and the other by $\mathbf{y}^{k+1}$ yields the update equation $\mathbf{y}^{k+1} = (I-AM^{-1})\mathbf{y}^k + \mathbf{b}$ which can also ...


5

Consider matrix a $2\times 2$ matrix $A$: $$ A=\left(\begin{array}{cc} 1 & 0\\ 2 & 1 \end{array}\right) $$ Singular values of $A$ are: $$ \sigma_1 = \sqrt{2}+1,\quad \sigma_2=\sqrt{2}-1 $$ resulting in the condition number (2-norm) $\kappa_2(A)=3+2\sqrt{2}$. If we consider a block-diagonal Jacobi preconditioner $J=I_2$ with block-size 1 (which ...


5

There are several things to consider in this experiment: Why Matlab sparse direct might be "so fast": (for your particular test) In 2D (of course, problem-dependent), your matrix $A$ arising after FEM discretization, after some reorderings might appear to be "close to banded" structure. The smaller is the bandwidth of $A$, the more efficient a sparse ...


5

Your argument is almost correct, but not quite. You correctly state that the use of one AMG step for $B^{-1}$ and the identity matrix for $S^{-1}$ leads to a constant number of iterations, but that is not the same as saying that it is optimal: If it takes me 1,000 iterations independent of the mesh size, then that's still a rather expensive method and I'd ...


4

Meijerink and van der Vorst showed that incomplete Cholesky dose not break down for $M$-matrices. As for finite element mass matrices, you will have to specify a basis to have any hope of making that claim. Given any basis $\Phi = [\phi_0 | \phi_1 | \phi_2 | \phi_3]$ such that the mass matrix $A = \Phi^T \Phi$ is SPD, we can construct a new basis $\hat \Phi ...


4

If you are only interested in the smallest eigenvalue, the conjugate gradient method applied to the matrix $L$ gives you a good approximation after a reasonably small number of steps, and you won't have to solve any linear systems. The details are in Y. Saad's book on iterative methods, but here a short summary: From the coefficients that are computed in ...


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