14

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


14

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


14

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

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

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

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

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

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


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

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

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

Performing $k$ steps of GMRES uses $O(n k^2)$ time and $O(n k)$ memory. In other words, the algorithm gets more and more expensive with each additional iteration. In theory, the algorithm terminates in $k \le n$ steps (ignoring round-off), thereby yielding a worst-case cubic $O(n^3)$ time and quadratic $O(n^2)$ memory figure. In practice, we manually ...


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


6

A few specific examples: Quasi-Newton methods avoid computing second derivatives by using an approximation to the Hessian which is updated at each iteration by a low rank (rank one or rank two) update. This makes factoring the Hessian (or equivalently keeping the inverse in product form) easy to do computationally. There are also limited memory Quasi-...


6

Given some $x_0$, you are looking for a bound on the error of the root for the equation: $ 0 = 0.7 \mathrm{sin}(x) - x + 5$ Note that the function is bounded above by $5.7 - x$ and below by $4.3 - x$ Each of these have a trivial root at $5.7$ and $4.3$ respectively. You know that any roots of your function then fall in $[ 4.3, 5.7]$. Therefore, for ...


6

I don't think there's enough information to give a heuristic like, "in general, it takes $k$ iterations for differential evolution to reach a global optimum". First, we don't know how many iterations deterministic global optimizers will take to converge to global $\varepsilon$-optimality in the general case. In limited special cases, we know that certain ...


6

Generally speaking, for many right-hand side (RHS) problems, a direct solver is a more feasible solution for several reasons: Major computations are performed during the factorization step (which is done only once for all RHS), and the solution find (for each RHS) is much cheaper. Direct solvers do not suffer from poor conditioning of the matrix or, in ...


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

The first thing you need to ask yourself: is your problem big enough that the overhead of MPI messaging is less than the work that you save. Your problem size is 10k which is small, but on the other hand you seem to have a dense matrix, so there is a good amount of work. Next: parallelizing a numerical method may change the mathematics. It looks like you're ...


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


6

It's a complicated question, which is why I've recorded a whole bunch of video lectures on the topic :-) Take a look at lectures 34 and following here: https://www.math.colostate.edu/~bangerth/videos.html


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 convergence of classical iterative solvers for linear systems is determined by the spectral radius of the iteration matrix, $\rho(\mathbf{G})$. For a general linear system, it is difficult to determine an optimal (or even good) SOR parameter due to the difficulty in determining the spectral radius of the iteration matrix. Below I have included many ...


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