47

Your question is a bit like asking for which screwdriver to choose depending on the drive (slot, Phillips, Torx, ...): Besides there being too many, the choice also depends on whether you want to just tighten one screw or assemble a whole set of library shelves. Nevertheless, in partial answer to your question, here are some of the issues you should keep in ...


16

Both of them are direct solver to solve linear systems (opposing to iterative solver). mldivide does perform the tests for $A$ in solving $Ax = b$. Please see Allan's answer in this thread for more information. Also see MATLAB's help on mldivide algorithm here. mldivide for square matrices: If A is symmetric and has real, positive diagonal elements, ...


15

Interesting that this question came yesterday, since I just finished an implementation yesterday that does this. My Background Just to start of, let me know that while my education background is from scientific computing, all work I have done since graduating, including my current Ph.D. work, has been in computational electromagnetics. So, I guess our ...


15

Have a look at Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods (Barrett et al.). You can find it here. Here's why I'm recommending this over other references: the "flowchart of iterative methods" in appendix D (last page) covers both linear solvers and preconditioners, it is short (100 pages or so), does not go into too ...


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

You can't beat an explicit formula. You can write down the formulas for the solution $x=A^{-1}b$ on a piece of paper. Let the compiler optimize things for you. Any other method will almost inevitably have if statements or for loops (e.g., for iterative methods) that will make your code slower than any straight line code.


13

When is a matrix ill conditioned? It depends on the accuracy of the solution you are looking for, as much as "beauty is in the eye of the beholder"... May be your question should better rephrased as are there cheap and robust condition number estimators based on the $LU$ factorization? Assuming you are interested in the real general (dense, non symmetric) ...


13

When you use ZGELSS to sovle this problem, you're using the truncated singular value decomposition to regularize this extremely ill-conditioned problem. it's important to understand that this library routine is not attempting to find a least squares solution to $Ax=b$, but rather it is attempting to balance finding a solution that minimizes $\| x \|$ ...


13

In general, there is no shortcut other than completely re-factoring the matrix. There have been a few similar questions on this SE that cover the topic in more depth than I can: Can diagonal plus fixed symmetric linear systems be solved in quadratic time after precomputation? LU Decom of PSD Matrix + Diagonal Matrix Perturbation of Cholesky decomposition ...


13

Ill conditioning is a feature of the system of the equations, not of the algorithm used to solve the system of equations. If your systems are that badly conditioned ($10^{15}$), then you can expect the solution to the system to be extremely sensitive to any perturbation of the problem data, even if the solution is done in extremely high precision (e.g. 500 ...


12

Your matrix $A$ isn't a circulant matrix- it's just Toeplitz. The method that you're trying to use fundamentally only works for circulant systems. Furthermore, your $a$ vector doesn't have the "-1" in it anywhere, so you clearly don't have sufficient information. A method that involves embedding the $n$ by $n$ Toeplitz matrix in a double-sized ...


12

I believe comparing an iterative method (multigrid) to a direct/exact method (Thomas) in terms of exact operation count isn't really meaningful. IIRC, Thomas operation count is $8N$ for any tridiagonal system. The only time I can imagine multigrid conceivably beating that is for a trivial case of having a linear solution, and even then the cost of evaluating ...


11

The short answer is that the Thomas algorithm will be faster than any iterative scheme for almost all cases. The exception would perhaps be applying a single iteration of a very simple iterative scheme such as Gauss-Seidel, but this is highly unlikely to give an acceptable solution. Also, this is ignoring parallel processing concerns. Multigrid is an ...


11

Practical experience shows that trying to get good initial iterates has little value. For example, in the context of solving partial differential equations, if you take the solution from one mesh, interpolate it onto a finer mesh, and use that as the starting guess for something like a CG iteration to solve the same problem on the finer mesh, then it turns ...


11

Defining the auxiliary variable $y=Bx$ yields the following algebraically equivalent expanded system, $$\underbrace{\begin{bmatrix} 0 & A \\ B & -I \end{bmatrix}}_{K} \underbrace{\begin{bmatrix} x \\ y \end{bmatrix}}_{u} = \underbrace{\begin{bmatrix} b \\ 0 \end{bmatrix}}_{f},$$ which you could solve with GMRES or another nonsymmetric Krylov method. ...


10

Without taking sides the discussion about whether to use direct or iterative solvers, I just want to add two points: There exist Krylov methods for systems with multiple right-hand sides (called block Krylov methods). As an added bonus, these often have faster convergence than standard Krylov methods since the Krylov space is built from a larger collection ...


10

Trivial answer for square $A$: use dgesvx which solves also for $A^T x = b$ when TRANS = 'T'. Please note that with BLAS or LAPACK you hardly have to transpose (swapping elements in memory) a matrix: most of the subroutines have a TRANS argument to accommodate for operation on the transpose matrix or on a matrix stored with a different memory layout. (...


10

This is called "structurally symmetric". It simplifies graph traversal, such as occurs when setting up aggregates in algebraic multigrid, but doesn't offer much structure to improve convergence rates. Note that all common PDE discretizations have this property so this is still a huge class of matrices including many instances for which no truly good ...


10

In some cases, (F)MG provides an algorithm with optimal properties. For instance, properly tuned FMG can solve some elliptic problems in a small number of "work units", where a work unit is defined to be the computational effort required to express the problem itself - in this case the operations to form the residual $b-Ax$ on the finest grid. This is such ...


10

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 matrix (e.g., symmetry), you could for example use GMRES. What you probably had in mind is the question of preconditioning, and that you can't use things such ...


9

There is typically a trade-off between the amount of work you put into constructing a good preconditioner for an iterative solver and the work you save by using a good preconditioner when actually solving the linear systems. In your case, the case is pretty clear: put as much work as you can into constructing a good preconditioner because you have to solve ...


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


9

Jed Brown has already pointed this out in the comments to the question, but there is really not very much you can do in usual double precision if your condition number is large: in most cases, you will likely not get a single digit of accuracy in your solution and, worse, you can't even tell because you can't accurately evaluate the residual corresponding to ...


9

There's no reason to append a row of 1's. You should just perform a rank-revealing QR factorization (like with routine SGEQP3) on $A^T$, and the last column of $Q$ should be in the nullspace. This has the added advantage that the relative magnitude of the last element on the diagonal of $R$ gives you some idea of how singular the solution is. Even better ...


9

Such a nested Krylov subspace method may work quite well in practice. It may be of interest for non symmetric linear systems for which restarted GMRES stagnates and unrestarted GMRES is too expensive or uses too much memory. Some literature: GMRESR: A family of nested GMRES methods, van der Vorst, Vuik Flexible inner-outer Krylov subspace methods, Simoncini,...


9

Since the matrix is so close to the identity, the following Neumann series will converge very rapidly: $$A^{-1} = \sum_{k=0}^\infty (I-A)^k$$ Depending on the accuracy required it might even be good enough to truncate after 2 terms: $$A^{-1} \approx I + (I - A) = 2I - A.$$ This might be slightly faster than a direct formula (as suggested in Wolfgang ...


9

Eigen 3 is a nice C++ template library some of whose routines are parallelized. c.f. Eigen documentation The parallelization is OMP only, so if you intend to parallelise using MPI (and OMP) it is probably not suitable for your purpose. The nice feature of Eigen is that you can swap in a high performance BLAS library (like MKL or OpenBLAS) for some routines ...


9

You want to minimize $\min \| Ax -y \|_{2}^{2} + x^{T}B^{T}Bx=\| Ax -y \|_{2}^{2} + \| Bx \|_{2}^{2}$ Recall that $\| u \|_{2}^{2} + \| v \|_{2}^{2}= \left\| \left[ \begin{array}{c} u \\ v \end{array} \right] \right\|_{2}^{2}$. Thus your problem can be written as $\min \| Hx - g \|_{2}^{2}$ where $H=\left[ \begin{array}{c} A \\ B \end{array} \...


8

The solution of an ill-conditioned system of equations with a matrix of norm 1 a random right hand side of norm 1 will have with high probability a norm of the order of the condition number. Thus computing a few such solutions will tell you what is going on.


8

For an unconstrained linear SPD problem, BFGS convergence matches unpreconditioned conjugate gradients. Differential operators appearing in structural mechanics problems are usually ill-conditioned with many large eigenvalues, making unpreconditioned iterative methods hopelessly slow. To deal with the troublesome spectrum, structural mechanics solvers ...


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