A lot of big sparse matrices

Hi everyone,

I need to get a sum of approx 10000 matrices, that are constructed by multiplication with sparse matrices.

using LinearAlgebra, SparseArrays
n = 101
p = 101
m = 25
#indK = sprand(m, n^2*p^2, 0.2)  EDIT on density
indK = sprand(m, n^2*p^2, 0.01) 
Sm = rand(m, m)
o = rand(n*p)
pom = zeros(n*p, n*p)
for j=1:n*p
     pom+= indK[:, (j-1)*n*p+1:j*n*p]'*Sm*indK[:, (j-1)*n*p+1:j*n*p]*o[j]
end

This is extremely slow, so I tranformed the line in for loop into

mul!(pom, indK[:, (j-1)*n*p+1:j*n*p]'*Sm, indK[:, (j-1)*n*p+1:j*n*p], o[j], 1)

but this is still too slow for me. This is part of another loop, indK does not change, but Sm and o do.
Any ideas on how to improve it or is this the best I can get?

Row-wise access is bad for CSC matrices. A lot of the numbers in the columns will likely be zero.

Note also that .2 is not sparse. Generally for sparse methods to be effective you need sparsity of .01 or so or lower. Also, using @views will make this much faster.

2 Likes

20% sparsity is probably too large to be worthwhile using sparse matrices. Does your real problem have more sparsity?

1 Like

The matrix in my original problem really is sparse, but its construction is too long to put it here.

Thanks for the note on sparsity, the matrix in my original problem is sparse, I will change the density.
Using @views as

@views pom += o[j]*indK[:, (j-1)*n*p+1:j*n*p]'*Sm*indK[:, (j-1)*n*p+1:j*n*p]

is way slower than the version with mul!, even when I convert the sparse matrix into dense.

Yeah, I was thinking about creating an array with matrices instead of one huge matrix, but it does not bring any improvement.

In that case, you lose the sparse data structure when you multiply by Sm, unless you use Sm = sparse(rand(m,m))… And pom is also not a sparse matrix.

3 Likes

Omg, you are right, how could I miss that. Together with taking indK as an array of sparse matrices, this solves my problem. Thanks a lot.