33

The important thing when choosing iterative solvers is the spectrum of the operator, see this paper. However, there are so many negative results, see this paper where no iterative solver wins for all problems and this paper in which they prove they can get any convergence curve for GMRES for any spectrum. Thus, it seems impossible to predict the behavior of ...


19

Introduce the vector $y:=-A^{-1}Gx$ and solve the large coupled system $Ay+Gx=0$, $G^Ty=-b$ for $(y,x)$ simultaneously, using an iterative method. If $A$ is symmetric (as seems likely though you don't state it explicitly) then the system is symmetric (but indefinite, though quasidefinite if $A$ is positive definite), which might help you to choose an ...


18

The choice between direct and iterative methods is dependent on goals and problem at hand. For Direct methods, we can note: The coefficient matrix of the linear system changes over the course of computation and may for sparse systems exhaust memory requirements and increase work effort due to fill-in Must complete to give useful results Factorization can ...


14

It depends a lot on the size of your matrix, in the large-scale case also on whether it is sparse, and on the accuracy you want to achieve. If your matrix is too large to allow a single factorization, and you need high accuracy, the Lanczsos algorithm is probably the fastest way. In the nonsymmetric case, the Arnoldi algorithm is needed, which is ...


14

You can just simulate the matrix-matrix product by forming the product of the two sparsity patterns -- i.e., you consider the sparsity pattern (that is stored in separate arrays in CSR format) as a matrix that contains either a zero or a one in each entry. Performing this simulated product only requires you to form the and operation on these zeros and ones ...


13

This is a well-studied problem in the field of sparse-direct solvers. I highly recommend reading Joseph Liu's overview of the multifrontal method in order to get a better idea of how reorderings and supernodes effect fill-in and solution time. Nested dissection is an extremely common way to generate the reordering, and essentially consists of recursive ...


13

I actually wrote the original code in Matlab for A*B, both A and B sparse. Pre-allocation of space for the result was indeed the interesting part. We observed what Godric points out -- that knowing the number of nonzeros in AB is as costly as computing AB. We did the initial implementaion of sparse Matlab around 1990, before the Edith Cohen paper that ...


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

BLAS1-operations, BLAS2-operations, and sparse-operations share the same curse of low arithmetic intensity, that they perform $O(1)$ flops for each memory read (contrast this to a BLAS3-operation like gemm, which performs $O(N^3)$ flops over $O(N^2)$ reads and only becomes more and more arithmetic-intensive/compute-bound for large $N$). However, sparse ...


12

The degeneracy of some eigenvalues looks to me like the hallmark of the breakdown of the Lanczos algorithm. The Lanczos algorithm is one of the more commonly used methods to approximate the eigenvalues and eigenvectors of Hermitian matrices; it's what scipy.eigsh() uses, through a call to the ARPACK library. In exact arithmetic, the Lanczos algorithm ...


12

The LU factors of a sparse matrix are at least somewhat sparse. The $Q$ matrix in QR can also somewhat preserve sparsity, and is typically used when the matrix is very long and skinny. The SVD of a sparse matrix will almost always have fully dense $U$ and $V$ factors, so it destroys any reason to perform the computations treating the matrix sparsely.


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

It's a question on what you spend your time on. For most of us, we spend 3/4 of the time programming and 1/4 of the time waiting for results. (Your numbers may vary, but I think the number is not completely without merits.) So, if you have a design that allowed you to program twice as fast (3/4 of a time unit instead of 1.5 time units), then you can can take ...


11

The cost of sparse matrix-vector multiplication scales linearly with the number of nonzero entries, as each entry is multiplied once by some entry in the vector. The cost of sparse matrix-matrix multiplication is highly dependent on the structure of the nonzeros. For instance, consider squaring a sparse matrix $A$ which is of an arrowhead structure: $$ A ...


11

For eigenvalues, simply take $k$ largest or smallest eigenvalues of $T$. They are good approximations of $A$, provided that the number of Lanczos iterations is large compared to $k$. Things are a little trickier if we want eigenvectors as well. The simplest way is to multiply each eigenvector $\mathbf{u}_i$ of $T$ by $V$ to the left, where $V$ is, as you ...


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

In deal.II (http://www.dealii.org -- disclaimer: I'm one of the principal authors of that library), we do eliminate whole rows and columns, and it is not too expensive overall. The trick is to use the fact that the sparsity pattern is typically symmetric, so you know which rows you need to look into when eliminating a whole column. The better approach, in ...


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


11

For what it is worth, for random sparse matrices of size 10,000 by 10,000 vs. dense matrices of the same size, on my Xeon workstation using MATLAB and Intel MKL as the BLAS, the sparse matrix-vector multiply was faster for densities of 15% or less. At 67% (as proposed by another answer), the dense matrix-vector multiplication was about three time faster.


10

All direct solvers supported by PETSc are available in Python under a common interface via petsc4py. Supported sparse direct solver packages include the PETSc native direct solvers, MUMPS, PasTiX, SuperLU, SuperLU_DIST, Umfpack, CHOLMOD, Spooles, LUSOL, MATLAB, and ESSL. See the MATSOLVER* man pages here.


10

PETSc uses BLAS for a few vector primitives, but these are generally limited by memory bandwidth and there isn't much variance in "optimization", so it tends not to make much performance difference. It also uses Lapack for some analysis such as Lanczos or Arnoldi estimates of eigenvalues and singular values, but these are generally not performance-sensitive....


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

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


9

A preconditioner, say M, is an approximation on the system matrix, say A that changes the problem into another problem with improved eigenvalue spectrum. A perfect preconditioner would be inverse of A i.e inv(M) = A. Unfortunately, this inverse is normally not avaiable, too complicated to compute, requires more space to store because of the fill-in's ...


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

Very few scientific software developers understand good principles of design, so I apologize if this answer is a bit long-winded. From a software engineering perspective, the goal of the scientific software developer is to design a solution that satisfies a set of constraints that are often conflicting. Here are some typical examples of these constraints, ...


8

SciPy supports sparse linear algebra via scipy.sparse.linalg (see the SciPy Documentation). SciPy supports the sparse direct solver packages SuperLU and UMFPACK.


8

I completely concur with the answers already given. I wanted to add that all iterative methods require some sort of initial guess. The quality of this initial guess can often affect the convergence rate of the method you choose. Methods like Jacobi, Gauss Seidel, and Successive Over Relaxation all work to iteratively "smoothen out" as much error as ...


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