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

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

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

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

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


7

The simplest/fastest way to solve ill-conditioned problems is to increase precision of computations (by brute force). Another (yet not always possible) way is to re-formulate your problem. You might need to use quadruple precision (34 decimal digits). Even though 20 digits will be lost in a course (because of condition number) you will still get 14 correct ...


7

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


7

One option here would be to form the normal equations $A^{T}Ax=A^{T}b$ and solve them by Cholesky factorization of the resulting $n$ by $n$ matrix. This squares the condition number of the problem which could potentially be a significant problem. Forming $B=A^{T}A$ doesn't require more than $O(n^2)$ memory, assuming that you can access the rows of $A$ ...


7

It's possible that your performance is limited by the memory bandwidth of your system. It's not at all uncommon to have a situation where two or three of your cores can use all of the available memory bandwidth and make it impossible to achieve higher parallel speedups because there isn't enough bandwidth to allow all of the cores to work at full speed. ...


7

There are a few ways you can do this. Even though RK4 fails unless stepsizes are really small, one thing that can happen a lot is that the model can start out more stiff than it really is in the full timespan. Thus I would still try something adaptive like MATLAB's ode45 or Julia's Tsit5() before ruling out non-stiff solvers. If that fails, then I would try ...


7

Several points I want to mention (with an encouragement to other CompSci users that are more familiar with Java specifics to give additional, more Java related answers): The solution of a system of linear equations and inversion of the matrix are two very different things. You almost never should explicitly invert the matrix. One should use one form of the ...


7

Using an Eigen matrix type where the number of rows and columns is encoded into the type at compile time gives you an edge over LAPACK, where the matrix size is known only at runtime. This extra information allows the compiler to do full or partial loop unrolling, eliminating lots of branch instructions. If you're looking at using an existing library rather ...


7

For a matrix that small, you're probably not going to do better than using dense methods. I wrote up a quick test in C++ for an 18x18 matrix with your sparse structure and randomly generated values and compared a sparse factorization (via Intel MKL's PARDISO) and the same dense factorization via Intel MKL's DGETRF. For the dense factorizations, any memory ...


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