MATLAB: Squaring a sparse upper-triangular matrix results in more non-zero elements than its original self?

MATLABsparse matrix

S=spconvert(dlmread('data.csv',','));
whos S: 400,000×400,000 double sparse 1.5GB
nnz(S): 100,000,000
numel(S): 1.6e11
S_2 = S*S;
whos S_2:400,000×400,000 double sparse 112GB
nnz(S_2) = 7e9 ????
numel(S_2) = 1.6e11
The original matrix S is sparse and strictly upper triangular, and the number of non-zero elements should therefore reduce in number every time it's multiplied by itself. The data for S is imported from CSV that has integer indices and values to 6d.p. Is this likely to be some floating point precision issue multiplying 0 to 6d.p by itself? I can't replicate the problem with small matrices, must be something in my data. I know it's big, but running it's no prob on a 512GB ram machine.

Best Answer

Why do you think the matrix should be more sparse after squaring?
Exactly the opposite should be the case, due to fill-in, at least if your matrix is truly sparse.
A = sprand(1000,1000,.1);
A = triu(A,1);
nnz(A)
ans =
       47715
nnz(A*A)
ans =
      401352
So I created A to be sparse, and strictly upper triangular, with zeros on the diagonal.
Perhaps you are thinking that since squaring such a strictly upper triangular matrix will now push the zeros one diagonal further up, So there will now be more zeros. But that is a small, even relatively tiny tradeoff compared to the fill-in that will happen.
So in the example above, I started with a randomly sparse strictly upper triangular matrix. The upper triangle was 10% non-zero. Square that, and the result will be almost completely filled in. One extra zero diagonal appears. But that is small potatoes.
Therefore unless your matrix was already virtually fully filled in, squaring it will result almost always in a more highly filled in matrix.
Only at the very end of the spectrum will there be a reduction in non-zeros. For example, a completely full upper triangle will indeed see a small reduction, as expected since the zeros propagate upwards.
nnz(triu(ones(1000),1))
ans =
      499500
nnz(triu(ones(1000),1)^2)
ans =
      498501
The difference lies in whether any fill-in will dominate the loss of a diagonal to zeros.
Related Question