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

S=spconvert(dlmread('data.csv',','));

whos S: 400,000x400,000 double sparse 1.5GB

nnz(S): 100,000,000

numel(S): 1.6e11

S_2 = S*S;

whos S_2:400,000x400,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.

Answer 1

[John D'Errico]

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.

[matt labber]

 Perhaps you are thinking that since squaring such a strictly upper triangular matrix will now push the zeros one diagonal further up

yes that's absolutely right! i'd been working with full upper triangular matrices all day, where the elements were pushed one away from the diagonal hence had fewer elements. user idiocy - please disregard.

Answer 2

[Jan]

If your input matrix is strictly upper triangular (all elements on and below the main diagonal are zero), squaring it will set the elements in the diagonal above the main diagonal to 0. Is it guaranteed that this decreases the number of zeros? Maybe there are many other zeros above the main diagonal which vanish at first, although the matrix is nil-potent (there is a positive k such that x^k is zero).

If (and only if) it is guaranteed, that squaring reduce the number of zeros, you have proven, that the input matrix is not strictly upper triangular. Check this again:

isequal(S, S-tril(S))

By the way: CSV files are text files. This is useful only if they are read or edited by human. for 100'000'000 non-zero elements and even more zeros this is a really bad format for storing data. The binary form or even better storing the sparse array in a MAT file would be much better.

[John D'Errico]

Actually, there are strictly upper triangular matrices for which squaring will reduce the number of non-zeros.

Squaring a strictly upper triangular matrix will migrate the zero diagonals upwards.

nnz(triu(ones(1000),1))
ans =
      499500

nnz(triu(ones(1000),1)^2)
ans =
      498501

It is clear that this matrix is strictly upper triangular. Of course, the number of non-zeros in the first matrix is easy to count:

1000*(1000-1)/2
ans =
      499500

After squaring, we can prove the number of non-zeros will be:

999*(999-1)/2
ans =
      498501

The reduction is small, but it is easy to prove the reduction.