MATLAB: Performance of matlab direct linear solver with respect to sparsity
backslash operatorsparse direct linear solver
Hi,
I'm using Matlab's direct linear solver (backslash operator) to solver a linear system generated from a mixed finite element discretization. Matlab says the matrix is sparse, when checked for sparsity (issparse()), however I know that the stiffness matrix is not that sparse compared to a conventianal finite element discretization. Now I have two questions,
I would like to know the effect of sparsity on the linear solver, lets say if the sparsity is increased, will the solver have the same performance, my benchmark for sparsity would be the sparsity of a stiffness matrix generated from a displacement based finite element.
I also see that the solver performes much better when I use full matrix (initiallise the stiffness matrix as zeros) compared to sparse. It looks like mutli-threading is fully exploited when full matrices are used (correct me if I'm wrong).
My stiffness matrix is not that big, close to 1e5.
Cheers,
ben
Best Answer
For sparse matrices, the performance of backslash depends a lot on the specific structure of that sparse matrices (where the nonzeros are placed).
In general, if a dense matrices was stored in sparse format, the backslash command would be slower for that than if it was in stored as a full matrix. This is because the sparse format is trying to exploit sparsity, but the matrix isn't actually sparse. In fact, a matrix needs to have significantly more zeros than nonzeros for a sparse solver to be faster than just calling the dense version.
So if you have tried this out and found that for your matrices, dense is faster than sparse, it's a good idea to stick to that. It's true that sparse algorithms are harder to thread than dense ones, so that can definitely be a factor - some of the sparse solvers do exploit threading, but in a less effective way than the dense ones.
Just to reiterate: The distribution of the nonzeros within a sparse matrix has a large impact on the performance of its solver, just knowing the size and the number of nonzeros isn't enough to predict the speed. The first step in a linear solver is to factorize the matrix (for example, compute triangular matrices L and U with A = L*U). For a sparse matrix, the L and U matrix might have many more nonzeros than the original A - this depends on the sparsity structure of A, and it's why that has a huge impact on performance.
The computational bottleneck in this code comes from converting an existing dense matrix into a sparse matrix. If the matrix has only a few non-zero entries, sparse matrix multiplication will indeed be much more efficient. However, there may be a large overhead cost in constructing the sparse matrix from the dense matrix that exists. The cost of the conversion process may therefore negate the savings gained in the multiplication.The best option is to construct the matrix A to be sparse from the beginning, rather than building it as a dense matrix and converting it:
>> A = sparse(i, j, s, m, n);
This will provide the performance gain of the sparse multiplication without the overhead cost of converting a dense matrix into a sparse matrix. However, if it s not possible to do so, then it could prove more efficient to simply perform a dense matrix multiplication.
You can define anything you want to be sparse if you so desire. So sparse(ones(1000)) will produce a sparse matrix. ;-)
But seriously, this matrix really is reasonably sparse. Just read what is shown on the page you direct to!
We see that it is shown to be an 1899x1899 square matrix. There are 20296 non-zeros in the matrix.
20296/1899
ans =
10.688
So on average, roughly 10.7 non-zeros per row of a matrix, where each row has 1899 elements.
20296/1899/1899
ans =
0.0056281
Total density of around 0.6% non-zeros. Is that matrix sparse? Well, yes, it is. Not massively so, but it is sparse. Since I don't have that matrix, I'll do a little memory comparison on a random one with the same density.
As = sprand(1899,1899,0.0056);
Af = full(As);
whos As Af
Name SizeBytesClassAttributes
Af 1899x189928849608double
As 1899x1899337568doublesparse
So we see that storing the matrix as sparse gives me a decrese in memory of almost a factor of 100-1. (Sparse matrix storage also requires we store the location of those elements, so it is not quite as efficient as we might like.)
b = rand(1899,1);
timeit(@() Af*b)
ans =
0.0023482
timeit(@() As*b)
ans =
6.1858e-05
And working with the matrix with a matrix multiply does give a significant savings.
Much of the time, when you work with a sparse matrix, you will perform a decomposition of some sort. Odds are that will not produce a speed increase here, because the matrix itself is not a nicely tightly banded matrix. I think much of the time, people think of a sparse matrix in the form of something like a tridiagonal matrix. A matrix factorization of a tridiagonal matrix will produce no fill-in at all. On the matrix you show, if we tried to do a Cholesky or LU factoriation, for example, the result will be a virtually completely full triangular matrix.
tic,[L,U] = lu(Af);toc
Elapsed time is 0.205769 seconds.
tic,[L,U] = lu(As);toc
Elapsed time is 3.032499 seconds.
As you see, not all computations on such a matrix will see a gain in speed. Here, we see the result ends up as a nearly full triangular matrix for the factors. So treating that matrix as sparse, IF you intend to do a decomposition of the matrix will be a time loss, because the overhead from the sparse algorithms is too costly.
nnz(L)
ans =
1177842
nnz(U)
ans =
1208969
numel(L)
ans =
3606201
That is to be expected on this matrix. It does not mean the original was not sparse. It is indeed sparse. Just perhaps not as massively sparse as you may want a matrix to be. Sometimes sparse must be measured by the eye of the beholder, in proper context.
Best Answer