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

8 次查看(过去 30 天)
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?

回答(1 个)

Steven Lord
Steven Lord 2023-5-10
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.
  3 个评论
Steven Lord
Steven Lord 2023-5-10
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
cTroels 2023-5-10
编辑:cTroels 2023-5-10
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.

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Sparse Matrices 的更多信息

产品


版本

R2023a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by