14

The most obvious thing you can do is to precompute [L,U] = lu(A) ~ O(n^3) Then you just compute x = U \ (L \ b) ~ O(2 n^2) This would reduce the cost enormously and make it faster. Accuracy would be the same.


14

There are two choices which, if you use reasonable data structures, can solve the problem with $n>10^6$ on a laptop and $n \approx 10^{12}$ on a supercomputer. Note that for efficiency, you should use multigrid to solve with $\Delta$. The cost in either case will be a small factor more expensive than just solving with $\Delta$. The two approaches are ...


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

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

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


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

We did some extensive computer labs in our scientific computing courses about this topic. For the "small" calculations we did there, Matlab's backslash operator was always faster than anything else, even after we had optimized our code as much as possible and re-orded all matrices beforehand (for instance with Reverse Cuthill McKee ordering for sparse ...


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

Suppose $A$ is a $n\times n$ dense matrix and you have to solve $Ax_i = b_i$, $i=1\dots m$. If $m$ is big enough then there is nothing wrong in V = inv(A); ... x = V*b; Flops are $O(n^3)$ for inv(A) and $O(n^2)$ for V*b, therefore in order to determine the break-even value for $m$ some experimentation is needed... >> n = 5000; >> A = randn(n,n)...


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


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

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


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


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