47

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


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

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


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

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


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

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


12

No, positive definiteness (and symmetry) are only precondition to using the Conjugate Gradient method. But there are plenty of other iterative methods such as MinRes and GMRES that can be used for indefinite and non-symmetric matrices.


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


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

Iterative Krylov-subspace solvers generally only require matrix-vector products and don't care whether or where there are zeros in the matrix. In your case, unless you have other information about the matrix (e.g., symmetry), you could for example use GMRES. What you probably had in mind is the question of preconditioning, and that you can't use things such ...


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


9

Yes you can do this, and it will converge in one iteration regardless of the starting value. This is because each step of Newton's method involves solving a linear system with the Jacobian of the nonlinear function. In this case the Jacobian just equals $A$. In other words: this is a little circular because it requires you to solve the system $Ax=b$ in the ...


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

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

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

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


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


7

If $q:=|f'(x^*)|<1$, where $x^*$ is the solution, the fixed point iteration you talk about is locally linearly convergent with convergence rate $q$. Thus if $q$ is small or zero, the method is competitive with Newton's method. Far away from the solution, convergence is difficult to predict in the absence of global information (such as a Lipschitz ...


7

Of course. The common Krylov space methods build the solution of $Ax=b$ as an approximation first in the space $\{r\}$ (where $r=b-Ax_0$ with the initial guess $x_0$), then $\{r,Ar\}$, then $\{r,Ar,A^2r\}$, etc. So the space in which the solution is sought depends on the right hand side vector. If your solution happens to be, say, in $\{r,Ar,A^2r\}$, then ...


7

Newton's Method can be interpreted as a type of fixed point iteration of the form $$x_{n+1}=F(x_n)$$ where $$F(x)=x-\frac{f(x)}{f'(x)}$$ Assume we know the root $r$ such that $F(r)=r$ (i.e. $r$ the answer to part 1 of your question). According to fixed point theory, if $F$ is a contractive mapping on a closed set $[a,b]$ and $F'$ exists and is ...


7

If your matrix is symmetric, positive definite, the CG method may converge slowly, but it converges for $n\to\infty$. The only reason it does not converge on a computer are round-off errors, in particular if the condition number of the matrix, the quotient of largest and smallest eigenvalue is large. Experience is, that with double precision arithmetic, ...


7

You can convert $\mathbf{b}-\mathbf{x}$ into polar coordinates, and do the dot product in this system. This changes $((\mathbf{b}-\mathbf{x})\cdot\nabla)^n\frac{1}{r}$ to $$\left((\mathbf{b}-\mathbf{x})_r\frac{\partial}{\partial r} + (\mathbf{b}-\mathbf{x})_{\theta}\frac{1}{r}\frac{\partial}{\partial \theta}\right)^n\frac{1}{r}$$ Here, I am using the ...


7

I've implemented this recently, basically it counts how many times each specific colour borders another colour to make up a frequency table. To generate an image, a random colour and position are selected and the rest of the image is built up from there. Results aren't very coherent, but they match the colour palette of the original. In order to make sure ...


7

These concepts are related. Let $A = M - N$ and consider the iteration $$ M x_{k+1} = N x_k + b. $$ We can write this as a mapping $\Phi : \mathbb{R}^n \to \mathbb{R}^n$ defined by $$ \Phi(y) = M^{-1}(N y + b), $$ and so $x_{k+1} = \Phi(x_k)$. $\Phi$ is called a contraction mapping if there exists a constant $0 \leq L < 1$ (called the Lipschitz constant, ...


7

We can't use the criterion you show last in practice because it requires us to know what $\kappa(A)$ is. But computing the condition number is, in general, more expensive than solving a linear system. As a consequence, the criterion you show is not practical. That only leaves us with variants of the criterion $$ \frac{\|r_k\|}{\|r_0\|}=\frac{\|b-Ax_k\|}{\|b-...


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