1
$\begingroup$

I’ve been facing a problem of solving linear systems Ax=b arising from discretized PDEs (Stokes equations in particular). Nively, it seems that solving Ax=b should not take much more time than simply multiplying matrix A by some vector v as long as A is sparse. More specifically, when discretizing a linear PDE, one end up with a sparce matrix A whose i:th row contains as many non-zero elements as there are nearest neighbors of node i within the spatial grid, and this is typically a small number. Thus, it seems that solving such a system should require O(n) time, where n is the size of the (square) matrix A. However, experience indicates that solving a linear system is way more costly than performing matrix multiplication, even in the case when A is sparse. In practice, I use MUMPS solver (that uses LU-factorization) implemented in PETSc library. I find that explicit solution becomes essentially unfeasible when the size of the matrix exceeds 100,000 or so. On the other hand, performing simple arithmetic manipulations with 100,000 floats seems like a trivial task. I come to the same conclusion experimenting with sparse matrices in MATLAB. For example:

A = sprand(100000,100000,0.0001) + speye(100000,100000);

v = ones(100000,1);

Now, multiplying A by v takes 0.02 seconds. However, solving A\v takes a much larger amount of time, even though A has roughly 10 elements on each row.

Is there a simple reason why solving a sparse system appears much harder than performing matrix-vector multiplication? Is there some solution to this?

$\endgroup$
4
  • $\begingroup$ Solving a linear system (even if it is sparse) costs much more time than a simple matrix-vector multiplication. In case you are using a direct solver, you have to perform the factorization as well as the forward and backward substitution. In case of you are using an iterative solver, you have multiple matrix-vector multiplication. I would recommend to read the theory behind solution of linear systems. $\endgroup$
    – vydesaster
    Commented Oct 29, 2021 at 21:21
  • $\begingroup$ Thanks vydesaster, you are basically acknowledging that I am not wrong in my statement. However, is there a concise reason for why this is so? I mean, my handwaving above seemingly suggests that the difference in speed should not be all that great... $\endgroup$
    – P. Trinli
    Commented Oct 29, 2021 at 21:43
  • $\begingroup$ I can't answer this in a satisfactory manner but I will try to tell you where you can go and find the answer. If you do something like a Gauss elimination you'll understand why. You'll start filling up your sparse matrix fairly quickly. You can also look at the LU decomposition algorithm. I am sure state of the art do something extra like reordering the matrix to minimize fill and something like that. The general rule is that if you are in 3D than you can only go up a few hundred thousand unknowns on your local machine. A bit higher in 2D because the matrix is more banded. $\endgroup$ Commented Oct 30, 2021 at 13:16
  • 1
    $\begingroup$ Maybe the main misunderstanding here is, that when you try to solve Ax=b, you are solving it for x. A and b are known. So there is no easy way to get x by a single matrix-vector multiplication. Instead you need to change the equation and bring A on the right side. This can be done theoretically by using the inverse of A. But this costs too much computation time. Thus you use a linear algebra solver. $\endgroup$
    – vydesaster
    Commented Oct 30, 2021 at 22:12

1 Answer 1

3
$\begingroup$

You can check Direct Methods for Sparse Matrices by Duff, Erisman and Reid (reviewed https://epubs.siam.org/doi/abs/10.1137/1031105) for the reason, but it won't be concise. For dense linear systems, analyzing the cost of multiplication and linsolve is relatively straightforward; you can even prove that the cost of linsolve is equal to matrix-matrix multiplication. These are nice results.

In case of sparse linear algebra, it is not that simple. You can prove that the cost of sparse matrix - DENSE vector multiplication is at O(nnz). But the moment you consider sparse matrix - sparse vector or sparse matrix - sparse matrix multiplications it gets tricky because now you have to consider cost of allocating new memory, cost of checking if you can do the operation $A_{ij}B_{jk}$ as $B_{jk}$ might be equal to zero, hence not stored at all, and a lot of other overhead.

Also, the standard way of doing sparse direct solves is to use LU factorization (or Cholesky or LDLT). Even if the matrix is sparse, its factors may not be. We try to guarantee this as much as possible using sparsity preserving permutations like COLAMD, or permutations along the way like Markowitz pivoting. Which is another cost.

Many other issues, like batch processing, cache misses, data movement... There are some special cases (like banded matrices with small bandwidths) on special (idealized) hardware where it is proven that the cost of solving a sparse linear system is O(nnz). That is it.

$\endgroup$
3
  • $\begingroup$ Thanks! If I understand your answer correctly, O(nnz) scaling for matrix vector multiplication implies that multiplying a matrix with e elements on every row by a dense vector takes ne time (since nz=e by definition). This is expected: one has to evaluate a linear combination of e numbers when looping every one of n rows. Now, by the same token, solving Ax=b by Gauss elimination should take ne^2 time which is not much more than ne. Why on Earth is this not the case?! $\endgroup$
    – P. Trinli
    Commented Oct 30, 2021 at 20:33
  • $\begingroup$ The cost analysis for the Gaussian elimination is not correct. At each step of the elimination, you may be, and this is very likely, creating new nonzeros. A classical example is the so-called arrowhead matrix, which is very sparse but its L and U factors are full. Which is why you need sparsity preserving (or fill-in reducing) permutations. Nevertheless, unlike the dense case, there is no straightforward relationship between SpMV and SpSolve operations. If you are into this, Tim Davis has his lectures on youtube. Give them a watch to see the complexity of sparse direct solvers $\endgroup$ Commented Oct 31, 2021 at 2:31
  • $\begingroup$ Thanks! - this rally clarifies things. $\endgroup$
    – P. Trinli
    Commented Oct 31, 2021 at 4:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.