(Very!) Slow elementwise division using large sparse matrices.

Hi!

I have written a pretty large code. A small part of the code takes the bulk of the solving time.

The small part of the code that takes a long time to solve is included here.

phi and x are N x 1 full vectors and H is a NxN sparse matrix. Sparse, both in the sense that it has a very small number of non-zero entries and it is the type of the matrix itself. N is approx. 50.000 and p is just a postive, real number.

sinxe = x.^p.*sin(phi);
cosxe = x.^p.*cos(phi);
Hssin = H*sinxe;
Hscos = H*cosxe;

Htsin = sinxe'.*H;
Htcos = cosxe'.*H;

HtDt = (Htcos./Hscos + Htsin.*(Hssin./Hscos.^2))./( Hssin.^2./Hscos.^2 + 1 )

The line with HtDt takes a long time to solve, even though im using sparse matrices and vectorized my inputs. Any ideas?

Answer 1

[Steven Lord]

The resulting matrix won't be very sparse. Anywhere your Hscos has a 0 (including the implicitly stored 0's) the result will contain a nonfinite, non-zero value.

Hscos = speye(5)

Hscos = 
   (1,1)        1
   (2,2)        1
   (3,3)        1
   (4,4)        1
   (5,5)        1

y = sparse(1)./Hscos

y = 
   (1,1)        1
   (2,1)      Inf
   (3,1)      Inf
   (4,1)      Inf
   (5,1)      Inf
   (1,2)      Inf
   (2,2)        1
   (3,2)      Inf
   (4,2)      Inf
   (5,2)      Inf
   (1,3)      Inf
   (2,3)      Inf
   (3,3)        1
   (4,3)      Inf
   (5,3)      Inf
   (1,4)      Inf
   (2,4)      Inf
   (3,4)      Inf
   (4,4)        1
   (5,4)      Inf
   (1,5)      Inf
   (2,5)      Inf
   (3,5)      Inf
   (4,5)      Inf
   (5,5)        1

fullY = full(y); % for comparison
whos Hscos y fullY

  Name       Size            Bytes  Class     Attributes

  Hscos      5x5               128  double    sparse    
  fullY      5x5               200  double              
  y          5x5               448  double    sparse    

You're consuming more memory storing y (every element of which is non-zero) as a sparse matrix than you would be if you stored it as a full matrix!

One potential approach would be to extract the non-zero elements from your sparse matrix Hscos, divide by the elements in the corresponding location in the other matrix, and treat anything with a 0 in Hscos as giving a 0 in the result. I don't know if that's the best approach for the underlying problem you're trying to solve, but if those 0's are actually important then you're going to need to operate on the full versions of your matrices.

[cTroels]

Thanks for the answer! Gave something to think about.. However, Hscos is not a sparse matrix, but a N x 1 vector of elements with the matrix-vector product of H and cosxe. Some of these elements may be very close to zero, but i find it unlikely that any of them are actually zero.

If N = 3, then the matrix HtDt would look like this. HtDt and H are banded matrices.

[Steven Lord]

Can you generate and show us what phi, x, and H look like for a small value of N (say N = 5 or N = 10)? And what is the order of magnitude for the elements of x and for p? I want to make sure you're not overflowing to Inf or underflowing to 0.

[cTroels]

Due to the current implemementation of the full code, i cannot (easily, atleast) compute it for N below 3600.

However, i can say that using spy(HtDt) that it shows a banded layout similar to H. Also, this size of HtDt is actually smaller than H.

x is between 1E-3 and 1, p is 3, phi is in the interval -pi to pi.

I know the entries in H is a finite real positive value in some interval i define myself and that it is banded along the diagonal.