15
votes
Accepted
Smallest eigenvalue without inverse
Compute the largest-magnitude eigenvalue $\lambda_{max}$ of $A$ (with, say, eigs('lm')).
Then compute the largest magnitude (negative) eigenvalue $\hat{\lambda}_{...
- 2,166
14
votes
Accepted
Quality of eigenvalue approximation in Lanczos method
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 ...
- 5,954
13
votes
Accepted
How does the QR algorithm applied to a real matrix returns complex eigenvalues?
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 (...
- 12.1k
13
votes
Accepted
When do not use preconditioners for sparse linear system of equations?
In my experience, you always need (or better use) some form of preconditioning. The type and complexity of the precondition would vary depending on the task though.
From Y. Saad, Iterative Methods for ...
- 8,592
12
votes
Accepted
Without positive definiteness, does an iterative solver work?
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 ...
- 53.7k
12
votes
Iteration counts of AMG solver changes in parallel
This is something that can happen in almost any numerical algorithm running in parallel.
It's important to know that floating-point addition is not associative due to round-off errors. Thus you can't ...
- 18.5k
11
votes
Accepted
Markov (Chain) image generators?
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 ...
- 226
10
votes
Accepted
How can a CG solver solve a non positive definite sparse matrix
I highly recommend the following read:
J.R. Shewchuk, "An Introduction to the Conjugate Gradient Method Without the Agonizing Pain"
In short, if the matrix is non-positive definite, there is no ...
- 8,592
10
votes
Accepted
Is there an iterative solver for dense matrices with possible zero diagonal entries?
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 ...
- 53.7k
9
votes
Accepted
Solve linear system with Newton-Raphson method
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 ...
- 3,263
9
votes
Accepted
Why MATLAB chooses the Householder in its built-in function gmres.m?
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 ...
- 2,069
8
votes
Accepted
How to solve the problem without using symbolic computation
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 ...
8
votes
Accepted
Efficient implementation of preconditioners for iterative solvers
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 ...
- 53.7k
8
votes
Accepted
What is the **contraction factor and convergence factor** of a iteration method?
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}(...
- 791
7
votes
Operation count for GMRES
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 ...
- 2,455
7
votes
Accepted
Why do not we choose the error solution norm as an iterative method's criterion?
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. ...
- 53.7k
6
votes
Accepted
Fastest way to solve a sparse unsymmetric system many times
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 ...
- 8,592
6
votes
For which problems Krylov subspace methods are preferred over multigrid methods?
The answer depends somewhat on the discretization; for example, some boundary integral discretizations result in very well-conditioned matrices, for which a Krylov solver works just fine without ...
- 3,102
6
votes
Accepted
Finite Elements: using preconditioned conjugate gradients with incomplete cholesky decomposition
While preconditioned CG with incomplete Cholesky (ICC) is reasonably straightforward to formulate mathematically, writing an efficient implementation is a non-trivial matter.
Here's some of the ...
- 2,166
6
votes
Accepted
Optimality of block-Jacobi preconditioner
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
$$
...
- 8,592
6
votes
Why my parallel code using MPI is much slower than the serial one?
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 ...
- 1,305
6
votes
Accepted
Iterative linear solver for "ugly" saddle point system
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 ...
- 2,069
6
votes
Accepted
Numerical Linear Algebra: When to use Direct methods versus iterative methods to solve a linear system - for PDEs in particular
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....
- 53.7k
6
votes
Why minimizing with respect to A-norm?
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 ...
- 10.6k
6
votes
Eigenvectors of Laplacian
They're on Wikipedia, for instance, in a page with the slightly unclear name of "Eigenvalues and eigenvectors of the second derivative".
- 10.6k
6
votes
Is there a way to generate a matrix-free decomposition for a matrix-free operator?
This is not efficient, primarily because (i) the factors of sparse matrices are, in general, not sparse themselves, and (ii) when computing matrix factorizations, you don't just need $A_{ij}$ but you ...
- 53.7k
5
votes
How to choose a method for solving linear equations
The Eigen Library Documentation also has a nice overview page with a lot of information about different matrix decompositions:
http://eigen.tuxfamily.org/dox/group__TopicLinearAlgebraDecompositions....
- 435
5
votes
MATLAB: code for restarted gmres
Technically, the template Matlab code provided at Netlib already calculates the obtained residual at every iteration. It's just not recorded in the way that this is an output of the subroutine.
After ...
- 8,592
5
votes
Accepted
Power Iteration over Rayleigh Quotient Iteration?
Power and Rayleigh iterations are fundamentally different beasts, because the former uses only matrix-vector products with $A$, while the latter uses matrix-vector products with $(A+zI)^{-1}$. This is ...
- 2,455
5
votes
Regarding impractical usage of direct solvers of linear systems
The "100,000 unknowns" rule-of-thumb applies to sparse matrices rather than dense ones. A naive direct solver, which doesn't take advantage of sparsity at all, could in principle have $O(n^3)$ ...
- 9,813
Only top scored, non community-wiki answers of a minimum length are eligible