New answers tagged iterative-method
6
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 ...
0
For Steel's answer. This is my matlab 2018b gmres.m: it seems that use Householder reflection. And your matlab 2018b gmres.m uses MGS.
Are you sure your matlab is 2018b? I am sure mine is matlab 2018b, firmly.
function [x,flag,relres,iter,resvec] = gmres(A,b,restart,tol,maxit,M1,M2,x,varargin)
%GMRES Generalized Minimum Residual Method.
% X = GMRES(A,B)...
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 ...
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 ...
1
Thanks for all your reply, I have done some experiments in matlab and confirmed that matlab indeed uses the sparse matrix-vector multiplication automatically optimally. For example,
clc;clear;
rng(0);
n=160;
A = gallery('poisson',n);% A is n^2 * n^2 and sparse matrix
x=rand(n^2,1);
tic;A*x;toc;
tic;full(A)*x;toc;
And the results are
Elapsed time is 0....
2
Multiplying matrix and a vector using sparse matrix format, such as CSR is not hard, and is a basic and common operation if you write a linear solvers library for sparse linear systems. Take SPARSEKIT of Yousef Saad as an example. In the code one of the first subroutines you will find is what he calls 'amux'. You can look at it below. As one can see it is ...
0
It does indeed work to linearize the forward mode this way, and will correspond to machine precision. The reverse or adjoint mode doesn't work this way despite my wishes to the contrary. The derivation I posted doesn't correspond exactly to machine precision, but the numerical tests I did showed that the magnitude of the error depends on the linear tolerance ...
2
Taking $x_0 = 0$, we have that $x_1 \in <M^{-1}b>$.
For the next iteration, we get $x_2 \in <M^{-1}b, M^{-1}AM^{-1}b>$.
For the next iteration, we get $x_3 \in <M^{-1}b, (M^{-1}A)^2M^{-1}b>$.
Continuing this argument, you will eventually find that $x_k \in \mathcal{K}_k(M^{-1}A,M^{-1}b)$.
This space can also be described as $Mx_k \in \...
3
The first version only converges if $\|I-A\|<1$.
The second version converges if $\|I-\omega A\|<1$, so the parameter $\omega$ allows you to use the iteration on a broader class of matrices.
The third version modifies $\omega$ at each step to minimize the $\|r_{k+1}\|_2$ with respect to $\omega_k$. The optimal value is given by $$\omega_k = \frac{...
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" ...
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-...
1
The catch here is that refining a mesh is easy/mechanical, but coarsening really isn't. So think about the problem from the other direction: mesh your problem as coarsely as possible, just enough to capture the (geometrical) features of the domain. This mesh is probably not fine enough to accurately represent the solution, but would be if you refined it ...
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 ...
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, ...
0
This statement seems a bit reductive for what is a rather large and involved problem. But multigrid, although it was developed for and is ideal for elliptic problems, is still just about the best we have for hyperbolic problems and sees widespread use (when people can manage to implement it). The comparison between Gauss-Seidel seems odd as they tend to ...
Top 50 recent answers are included
