45

Your question is a bit like asking for which screwdriver to choose depending on the drive (slot, Phillips, Torx, ...): Besides there being too many, the choice also depends on whether you want to just tighten one screw or assemble a whole set of library shelves. Nevertheless, in partial answer to your question, here are some of the issues you should keep in ...


25

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


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

I originally didn't want to give an answer because this deserves a very long treatment, and hopefully someone else will still give it. However, I can certainly give a very brief overview of the recommended approach: Perform a thorough literature search. If that fails, try every preconditioner that makes sense that you can get your hands on. MATLAB, PETSc, ...


15

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

Others have already commented on the issue of preconditioning what I will call "monolithic" matrices, i.e. for example the discretized form of a scalar equation such as the Laplace equation, the Helmholtz equation or, if you want to generalize it, the vector-valued elasticity equation. For these things, it is clear that multigrid (either algebraic or ...


13

This list is nowhere near complete, but hopefully the size of it will give a hint as to the scale of possible factors. I am assuming you are compiling the code from source on your platform of choice. Software Standard Library Performance Lin. Alg. Library Performance (if the software links to outside libraries) Compiler Choice Compiler Optimization ...


13

Compute the largest-magnitude eigenvalue $\lambda_{max}$ of $A$ (with, say, eigs('lm')). Then compute the largest magnitude (negative) eigenvalue $\hat{\lambda}_{max}$ of $M = A - \lambda_{max}I$ (again, through a standard call to eigs('lm')). Observe that $\hat{\lambda}_{max} + \lambda_{\max} = \lambda_{min}(A)$. The reason why this holds is explained here. ...


13

In a nutshell, the QR algorithm applied to a matrix $A$ is an iterative procedure that converges to the real Schur decomposition: a unitary matrix $Q$ and a matrix $R$ in block upper triangular form (see below) such that $A = QRQ^T$. It follows that the columns of $Q$ are the eigenvectors (which are the principal objects that are computed!) and that $R$ has ...


12

Jack has given a good procedure for finding a preconditioner. I will try an address the question, "What makes a good preconditioner?". The operational definition is: A Preconditioner M accelerates the iterative solution of $A x = b$, and $M^{-1}$ can be applied cheaply compared to $A^{-1}$. however this does not give us any insight into designing a ...


12

Damped Jacobi Suppose the matrix $A$ has diagonal $D$. If the spectrum of $D^{-1}A$ lies in the interval $[a,b]$ of the positive real axis, then the iteration matrix of Jacobi with damping factor $\omega$ $$B_\text{Jacobi} = I - \omega D^{-1} A$$ has spectrum in the range $[1 - \omega b,1 - \omega a]$, so minimizing the spectral radius with $$\omega_{\text{...


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

The derivation of the BFGS is more intuitive when one considers (strictly) convex cost functionals: However, some background information is necessary: Assume, one wants to minimize a convex functional $$ f(x) \to \min_{x\in \mathbb R^n}. $$ Say there is an approximate solution $x_k$. Then, one approximates the minimum of $f$ by the minimum of the ...


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

In exact arithmetic you shouldn't need to reorthogonalize regularly, but practically you do. Your u1 and u2 are close to (but not exactly) the true eigenvectors, so your initial deflation almost (but not entirely) removed the true eigenvectors from u3. The tiny components you left behind will be amplified by repeated multiplication by A, you will need to ...


11

In addition to Christian's answer, it's also worth noting that for linear convergence you have $e_{k+1} \le \lambda_1 e_k$ where you have $\lambda_1<1$ if the method converges. On the other hand, for quadratic convergence you have $e_{k+1} \le \lambda_2 e_k^2$ and the fact that a method converges does not necessarily imply that $\lambda_2$ must be smaller ...


11

You define a sequence in, say, $C^{\infty}(\Omega)\times C^{\infty}(\Omega)$ by $$ \frac{\text{d}^2u^k}{\text{d}x^2} + \frac{\text{d} v^{k-1}}{\text{d} x} =f\\ \frac{\text{d}^2v^k}{\text{d}x^2} + \frac{\text{d} u^{k-1}}{\text{d} x} =g\\ $$ (plus boundary conditions). It is clear that if this sequence converges, it will be a solution of your original set of ...


11

In general, there is no shortcut other than completely re-factoring the matrix. There have been a few similar questions on this SE that cover the topic in more depth than I can: Can diagonal plus fixed symmetric linear systems be solved in quadratic time after precomputation? LU Decom of PSD Matrix + Diagonal Matrix Perturbation of Cholesky decomposition ...


10

It is critical to know more about the structure. It matters whether the random entries are uniformly or normally distributed and whether there is a shift or not. If there is no structure at all, then you cannot asymptotically beat a direct solve. Some comments on your proposed approaches Incomplete LU is complete LU when applied to a dense matrix. You could ...


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

The convergence behavior you are seeing is actually expected. One of things that makes the Lanczos method so interesting is that it does a good job of simultaneously converging eigenvalues at both ends of the spectrum. I assume your expectation of converging only the largest eigenvalues is based on the fact that, as expected from the Power iteration ...


9

For a single rational equation in the complex domain, the basin of attraction is fractal, the compelement of a so-called Julia set. http://en.wikipedia.org/wiki/Julia_set . For theory with some nice online figures, see, e.g., http://mathlab.mathlab.sunysb.edu/~scott/Papers/Newton/Published.pdf http://hera.ugr.es/doi/15019160.pdf Even the ''globalized'' ...


9

In practice, yes. While $e_k$ is still large, the rate coefficient $\lambda$ will dominate the error rather than the q-rate. (Note that these are asymptotic rates, so the statements you linked to only hold for the limit as $k\to\infty$.) For example, for first order methods in optimization you often observe an initially fast decrease in error, which then ...


8

The conjugate gradient algorithm works for semidefinite problems and produces the minimal norm solution.


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

Youssef Saad's book Numerical Methods for Large Eigenvalue Problems, 2nd edition uses the norm of the residual vector to define convergence criteria. He defines the residual vector as follows on page 59: Given a matrix $\mathbf{A} \in \mathbb{C}^{n \times n}$, a putative eigenvalue $\widetilde{\lambda} \in \mathbb{C}$ and a putative eigenvector $\widetilde{\...


8

There is a theorem that syas that a black box algorithm is guaranteed to find the global minimum of an arbitrary smooth (i.e., twice continuously differentiable) function if and only if it samples points densely in the search space. Here dense is meant in the topological sense, i.e., it must sample points in arbirarily small neighborhoods of every point. ...


8

In theory, yes. In practice, rounding errors will usually result in (initially slow) convergence to $u_1$. At essentially the same cost one can run the Lanczos algorithm, which will have much faster convergence, and produce the three dominant eigenvalues unless two of these eigenvalues are essentially the same. For Lanczos, selective reorthogonalization is ...


8

The simple answer is that you would use inverse iteration (subspace or with deflation). This is basically the power method (repeatedly multiplying the matrix by a vector and normalizing, singling out the eigenvector corresponding to the eigenvalue of largest magnitude) applied to $A^{-1}$. Since you desire the $k$ closest to the origin, you need to use some ...


8

You can solve this numerically in Python without symbolic computation. from __future__ import print_function, division import numpy as np from numpy import exp from scipy.integrate import quad from scipy.optimize import root def f1(a1, a2, x): return exp(a1 * x + a2 * x * x * x) / (1 + x * x) def f2(a1, a2, x): return exp(a1 * x + a2 * x * x * x) *...


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