26

Yes, the condition number always matters in floating-point arithmetic, whether you choose to solve your system with an iterative or direct method. The relative accuracy of an approximate solution to $Ax = b$ obtained from LU factorization with pivoting is $O(\kappa(A) \cdot \varepsilon)$, where $\varepsilon$ is the smallest floating point number such that $1 ...


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

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


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

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


12

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


12

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


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


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

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

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

In the context of finite element methods (and, especially, symmetric problems) the most common direct solution method is Cholesky factorization (plus following substitutions). MATLAB uses Tim Davis' CHOLMOD package to compute Cholesky factorization whenever the heuristics of backslash operator encounter a symmetric positive definite matrix. In fact, Julia ...


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

One thing that CG has in its favor is that it's not minimizing the discrete $l^2$ norm for its residual polynomials (what GMRES does). It's minimizing a matrix-induced norm instead, and very often this matrix-induced norm ends up being very close to the energy norm for discretizations of physical problems, and frequently this is a much more reasonable norm ...


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


8

There is an extensive literature on numerical aspects of the implementation of the simplex method including methods for updating the factorization of the basis as well as techniques that maintain the inverse of the basis in product form. Replacing a column of the basis matrix is a "rank-one update" since it can be written in terms of adding a matrix to the ...


8

We should be more precise here. The simplest estimate that you can give is that $$ ||x^* - x|| = ||A^{-1} A (x^* - x)|| \le ||A^{-1}||\,||b - A x|| $$ so that if you terminate your iteration using the residual, you can be off by a factor of $||A^{-1}||$, for the relative residual by $$ \kappa(A) = ||A||\,||A^{-1}|| $$ so you have a simple estimate of how ...


8

Your two ideas make it much too complicated. If $X$ is the inverse of $A$, $$ AX=I, $$ and $x_i$ is the $i$-th column of $X$ and $e_i$ is the $i$-th column of the identity matrix $I$ ($e_i$ is a vector of all zeros except with $1$ in the $i$-th location), then the columns $x_i$ of the inverse are defined by $$ Ax_i = e_i. $$ All you need to do is solve $n$ ...


8

The matrix math in this paper is terribly hard to follow, but here's what equation (19) should look like $\left( \begin{array}{c|c|c|c} c_{11} \mathbf{I} - \mathbf{A}^T & c_{12} \mathbf{I} & c_{12} \mathbf{I} & \mathbf{0}^{3\times 3} \\ \hline c_{21} \mathbf{I} & c_{22} \mathbf{I} - \mathbf{A}^T & c_{23} \mathbf{I} & \mathbf{...


8

It is standard procedure for general case like that. You can enforce condition for electrode by Lagrange multipliers, $$ L(u,\lambda) = \int_\Omega \sigma u_{,i} u_{,i} \, \textrm{d}V - \int_{\partial\Omega_\sigma} u n_i j_i \textrm{d}S + \lambda\left(\int_E \sigma u_{,n}\, \textrm{d}S - I \right) $$ that for discretised problem is $$ \left[ \begin{array}{...


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