MATLAB: Problem of “sprandn”: the number of nonzero elements is inconsistent with density
sprandn
R = sprandn(m,n,density) is a random, m-by-n, sparse matrix with approximately density*m*n normally distributed nonzero entries (0 <= density <= 1).
However, if density is large, the number of nonzero entries it not consistent with density.
nnz(sprandn(100,10,0.3))
ans =
259
nnz(sprandn(100,10,0.7))
ans =
492
For density=0.7, the number of nonzero elements should be approximately 700, but the result is only 492.
This is an example, I have tried for many times but results are similar.
I apologize in advance if it is a trivial problem.
Best Answer
The approximation is only good when density << 1. I think SPRAND/SPRANDN simply create random draw of (density*m*n) pairs of (i,j) WITH replacement (perhaps for the reason of efficiency and memory).
When density specified is close to 1, there are a lot of pairs that collide in the same place and the effective density is less than what user specifies.
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.
A determinant is NEVER a good measure of whether a matrix is singular. And as matrices get large, the determinant becomes more than useless.
For example,
det(eye(1000)) == 1
however,
det(0.1*eye(1000)) == ???
When computed in floating point in MATLAB, the latter will yield zero, which I imagine you know not to be true. However, both matrices are equally non-singular. As easily, I could have modified the problem to yield inf by a different scale factor.
Answer? Don't use the determinant for anything of importance. It will only guide you to a poor conclusion.
Of course, this does not imply that your matrix has been created correctly or incorrectly. However, your test of that using the determinant is a useless one. Instead, you might look at the rank, or the eigenvalues. A valid stiffness matrix will be at least positive semi-definite. Don't forget that even for a singular matrix, eig can return negative eigenvalues on the order of -eps.
Best Answer