# Tag Info

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.

7

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

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

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

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

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

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

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

A simple approach: F = eye(7); % gives you the Kronecker delta part for mu = -3:3 for nu = -3:3 isum = 0; % temporary variable for the sum for j = 0:19 isum = isum + besselh(mu-nu,7.0*j); end F(4+mu,4+nu) = F(4+mu,4+nu) + (1i)^(nu-mu) * isum; end end For systems larger than 7x7, it might be wise to precompute the bessel function ...

3

LSQR uses only 4 N-vectors of MEMORY. It might be your choice if storage must be minimized. It is fairly good otherwise, too, but probably not the best one (unless your matrix is non-square; then the ones below do not work bot LSQR does). LSQR is CG for $A^*A$ but numerically more robust etc. 1DOTs/MatVec, 3AXPYs/MatVec (2 MatVec). https://web.stanford.edu/...

3

Since the matrix is considerably sparse and well-conditioned (if it is true), I suggest you can try to use Krylov subspace method, which can only use the information of matrix-vector product; such as GMRES(m), CGS, BiCGSTAB, and IDR(s). With aspect of memory and well-conditioned matrix, maybe BiCGSTAB, CGS and IDR(s) are preferable. They are many FORTRAN ...

3

It seems you are looking for an optimization (here, specifically a fit). You have some experimental data and you want to have a simple model explaining it. The simplest way is using the Ordinary Least Squares. The idea is to minimize the standard deviation or the error you make on the vertical axis. Another one would be the PCA where you try to find an ...

3

The eta factorization of the basis is a technique widely used in the simplex method for LP to handle rank one column updates. The same idea can be modified to handle rank one row updates, I'll assume that you start with an LU factorization of $A_{0}$ and can solve systems of equations of the form $A_{0}x=b$. We can also solve systems of the form $yA_{0}=c$...

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