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.


9

The MUMPS sparse direct solver can handle symmetric indefinite systems and is freely available (http://graal.ens-lyon.fr/MUMPS/). Ian Duff was one of the authors of both MUMPS and MA57 so the algorithms have many similarities. MUMPS was designed for distributed-memory parallel computers but it also works well on single-processor machines. If you link it ...


7

You may want to watch lecture 34 here: http://www.math.tamu.edu/~bangerth/videos.html


7

The general problem that direct solvers are suffering from is the fill-in phenomenon, meaning that the inverse of a sparse matrix may be dense. This leads to huge memory requirements if the structure of the matrix is not "suitable". There are attempts to work around these issues, and MATLAB's default lu-function employs a few of them. See http://...


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


7

You are effectively asking how to compute $y=(\log M )^{-1}b$, where $M=e^A$ is the given matrix. There are several methods for computing $f(M)b$ without forming $f(M)$, and they are reviewed here. One general method is to use Cauchy’s theorem, $$y=\dfrac{1}{2\pi i}\int_\Gamma f(z)(zI - M)^{-1}b\,dz,$$ with $f(x) = 1/\log(x)$. $\Gamma$ is a contour that ...


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

Prof. Bangerth's lecture has already covered most of it, but I'd recommend you look at this paper. The authors take three of the most common methods (GMRES, CGS and CGNE) and give some matrices for which one method will converge in O(1) operations and the other two converge in O(N). The upshot of this is that there is no uniformly best iterative method for ...


6

The condition number for the stiffness matrix of any method (finite elements, finite volumes, finite differences) applied to a second order differential operator always grows as ${\cal O}(h^{-2})$ where the exponent equals the order of the differential operator. As a consequence, if you just make the mesh fine enough, you can make the condition number ...


6

Your second term is linear in $x$, so it can be rewritten as $Cx$ where $C$ is a suitable $n\times n$ matrix. You just need to figure out what its entries $C_{ij}$ are: for this you need a little index manipulation, but this should be an easy exercise in linear algebra. Alternatively, you can use the (much better) notation $B(x\otimes v)$ for the second ...


6

Based on Federico's answer, I think you can easily figure out what are $C_{ij}$ entries by explicitly writing the second term ($Bxv$) as: $$(Bx)_{ij} = \sum_{\alpha=1}^{n} B_{ij\alpha} x_{\alpha}$$ $$(Bxv)_{i} = \sum_{\beta=1}^{n} (Bx)_{i\beta} v_{\beta}$$ $$(Bxv)_{i} = \sum_{\beta=1}^{n} \sum_{\alpha=1}^{n} B_{i\beta\alpha} x_{\alpha} v_{\beta} = \sum_{\...


6

As Federico has mentioned, you probably don't want to deprive yourself of the learning experience. I'll just give you a small nudge in the right direction. You will never be able to store $A$. You also won't be able to store $(A+uv^T)^{-1}$. However, you don't really need to. You can easily write down a formula for each of the entries in $A$. Instead of ...


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


5

While it is easy to reach a suboptimal solution for your problem it is usually much harder (and problem-dependent) to come up with an optimal/robust strategy. In this sense your question is quite broad and to give good advice it is crucial to know about the problem you want to solve. In the following I will assume you speak of sparse systems and you have ...


5

You've started with a singular linear system of equations $Ax=b$. As a practical matter, it's unlikely that $b$ lies exactly in the range of $A$, so at best you can find a least squares solution that minimizes $\min \| Ax - b \|_{2}$ Because the system is singular, the null space of $A$ is non-empty, and there will be an infinite number of solutions to ...


5

This kind of scaling is fairly common and sparse direct factorization methods are commonly used on matrices of up to a few hundred thousand rows and columns. By the time you get to $N=20$, you’ll have a million rows with a thousand nonzero entries per column. This is getting into the range where sparse direct factorization becomes impractical and you’ll ...


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


5

First of all, MATLAB's gmres assumes that the preconditioner you use is linear. This is important! Actually it is the main difference between FGMRES and GMRES. Right preconditioned GMRES and FGMRES are exactly the same if you use a linear preconditioner, however, FGMRES allows the use of non-linear preconditioners. What do I mean by a non-linear ...


4

The closed-form solution is available by projecting the problem into the space spanned by B. To see this, note that we have $$\min_{v\in\mathbb{R}^{n}}\|X-BD_{v}B^{T}\|_{F},$$ but if we introduce $\tilde{X}$ such that $X=B\tilde{X}B^{T}$, then the optimization is reduced $$\arg\min_{v\in\mathbb{R}^{n}}\|B\tilde{X}B^{T}-BD_{v}B^{T}\|_{F}\equiv\arg\min_{v\in\...


4

It appears to me that you are dealing with a sequence of linear systems $A_j x_j = b$, where $A_{j+1}$ is a low rank modification of $A_{j}$. In your case, I would investigate if the Krylov subspace $K = \text{range $V$}$ which you have built to solve one linear system, say, $A x = b$ is relevant for the solution of the next problem, i.e., $$(A + \Delta A) ...


4

Michael Saunders wrote a sparse LU package that can do rank-1 updates, LUSOL. You could try to use that, since you write that direct solvers are viable for your problem.


4

Welcome to the site. You are actually virtually finished, but may not have realised it yet. A tridiagonal linear system is another name for a matrix problem which only has non zero entries on the leading diagonal and the one above and below it, so lets try writing your problem like that. We want a form $$ \mathbf{A} (\mathbf{\Delta W}^m) = \mathbf{b},$$ ...


4

In my opinion, you need a really good excuse to take a symmetric positive definite system and then turn it into an indefinite system. The first thing I would try is the conjugate gradient method with your original $R^{T}R+D=R^{T}b$ system, except perform 3 matrix-vector products each iteration. Treat $R^{T}R+D$ as an operator instead of explicitly forming ...


4

You would need to linearize the problem. I prefer to do it before discretization but it's possible to do also after discretization. (I'm a bit skeptical of linearization after discretization because I have never looked into the details. In general, discretization and linearization steps do not commute.) In the following I assume that the equation is actually ...


4

Most people use (algebraic as well as geometric) multigrid as preconditioners these days. It's an empirical observation that that leads to faster convergence in terms of iterations, given that in a typical CG iteration, applying a multigrid preconditioner takes up the vast majority of the effort compared to all other operations. In other words, you get the ...


4

I'd recommend a BLAS1/2-style implementation of LU decomposition with partial pivoting, along with forward/backward triangular solves. If you only have single digit numbers of unknowns, it should be fine to run dense/O(N^3) algorithms even on structured systems (ie banded/sparse/?). All these algorithms can be written in terms of "axpy" operations, ...


3

Ultimately, it depends on the sparsity of $A$ and $B$ and the symmetry of the resulting hadamard product. The hadamard product will aggregate the sparsity structure of $A$ and $B$. So if one or the other is sparse, the product is also sparse (and may be more so if $A$ and $B$ have difference sparsity structures). So any sparse direct solver would be ...


3

If you have access to the MATLAB optimization toolbox then this can easily be done using the quadprog() function. You'd start by writing the objective in quadratic form as $ \| Ax - b \|_{2}^{2} = x^{T}(A^{T}A)x-2(A^{T}b)^{T}x+b^{T}b$ then multiply to get $P=A^{T}A$ and $q=-2A^{T}b$. Then your objective is $f(x)=x^{T}Px+q^{T}x+b^{T}b$ and ready to ...


3

Direct vs. iterative is certainly one of the key questions in linear solvers for sparse matrices as you have already observed. There are many misconceptions and lots of misinformation on this topic. So I encourage you to approach the issue from this point of view rather than looking for "5 methods." Here is a nice introduction from one of the major ...


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