# Elementary question on numerical linear algebra

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?

• 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. Commented Oct 29, 2021 at 21:21
• 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... Commented Oct 29, 2021 at 21:43
• 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. Commented Oct 30, 2021 at 13:16
• 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. Commented Oct 30, 2021 at 22:12

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.