Can I reduce computation time for the matrix computation (sparse matrices)

Hello, in the code below, I have to do matrix multiplication thousands of time. The time it takes on my PC is about 70 secs to run this code. I want to reduce the time of execution. In the code below two variables, PRODOTTO3 and COEFF are sparse matrices and sparse vector respectively. I have tried making them sparse matrix but the time of execution increase to more than 200 secs. Does anybody have better idea to reduce the execution time?

From the time profiler I know that the 99.8% of the execution time is spent on this line "COEFF_S(i,j) = sum(COEFF.*(PRODOTTO3(:,j,i)))"

I really want to reduce the execution time.

THIS IS THE CODE WITHOUT SPARSE MATRIX

NU =25; %No. of Uncertain input Parameter
p=2 ; %Keep it same, 
npolinomi=factorial(NU+p)/(factorial(NU)*factorial(p)); % No. of PCE terms (Polynomials)  Formula is (N+P)! / (N!P!)
PRODOTTO3 = (double((randn(round(npolinomi),round(npolinomi),round(npolinomi)))>0.7));
COEFF = zeros(round(npolinomi),1) ; 
COEFF(1:3) = 1 ; 
COEFF_S = zeros(size(COEFF,1),size(COEFF,1)) ; 

tic ; 
for l = 1:200 %freq.
for i  = 1 : size(COEFF,1)
  for j=1:size(COEFF,1)
  COEFF_S(i,j) = sum(COEFF.*(PRODOTTO3(:,j,i)))  ;
end 
end
end%freq.

for i= 1 : 250 %freq
  for j = 1 : K
      COEFF_S = COEFF_S*COEFF_S ; 
  end
end
end %freq
toc

THIS IS CODE WITH SPARSE MATRIX

if true 
NU =25; %No. of Uncertain input Parameter
p=2 ; %Keep it same, 
npolinomi=factorial(NU+p)/(factorial(NU)*factorial(p)); % No. of PCE terms (Polynomials)  Formula is (N+P)! / (N!P!)
PRODOTTO3 = ndSparse(double((randn(round(npolinomi),round(npolinomi),round(npolinomi)))>0.7));
COEFF = zeros(round(npolinomi),1) ; 
COEFF(1:3) = 1 ; 
COEFF=sparse(COEFF) ; 
COEFF_S = zeros(size(COEFF,1),size(COEFF,1)) ; 

tic ; 
for l = 1:200 %freq.
for i  = 1 : size(COEFF,1)
  for j=1:size(COEFF,1)
  *COEFF_S(i,j) = sum(COEFF.*(PRODOTTO3(:,j,i)))  ;*
end 
end
end%freq.

for i= 1 : 250 %freq
  for j = 1 : K
      COEFF_S = COEFF_S*COEFF_S ; 
  end
end
end %freq
toc
end

[Jan]

What is K?

[Abi Waqas]

For you comment on this piece of code, originally it looks like this for K =1, it will be A(:,:,i), then K=2, it would be B(:,:,i), so forth and so on. In short for every value of K, two matrices in multiplication are changing. How can I reduce the execution time of this?

if true

 for i= 1 : 250
    for j = 1 : K
        COEFF_S(:,:,i) = COEFF_S(:,:,i)*A(:,:,i) ; 
    end
end
  end

Answer 1

[Matt J]

This double loop

    for i  = 1 : size(COEFF,1)
      for j=1:size(COEFF,1)
       COEFF_S(i,j) = sum(COEFF.*(PRODOTTO3(:,j,i)))  ;
      end
    end

should be replaced by

COEFF_S=reshape( COEFF.'*PRODOTTO3(:,:) , size(COEFF,1), []).';

This loop

    for j = 1 : K
      COEFF_S = COEFF_S*COEFF_S ; 
   end

should be replaced by

     COEFF_S=COEFF_S^K;

[Abi Waqas]

Dear Matt J Thank you for the answer. Second suggestion in not valid in my case because matrices are chaging everytime. I didn't show this for simplicity. I can't skip this loop for technical reason.

For the first suggestion I will definitely try and get back to you.

What about the exploitation of sparse matrix? Majority of the element of variable PRODOTTO3 and COEFF are zeros and I want to exploit that information to reduce the time of execution.

[Abi Waqas]

Dear Matt J

Thank you so much for the suggestion that definitely reduces the time of execution. But look at the code below, when I make PRODOTTO3 and COEFF sparse the execution time increases. Why is this happening? I really want to exploit the information both variables are sparse.

 if true
 PRODOTTO3 = ndSparse(double(351,351,351)>0.7));

 COEFF = zeros(351,1) ; 
 COEFF(1:3) = 1 ; 
 COEFF=sparse(COEFF) ; 
 COEFF_S = sparse(size(COEFF,1),size(COEFF,1)) ; 
 tic ; 
 for l = 1:200 %freq.
 COEFF_S=reshape( COEFF.'*(PRODOTTO3(:,:)) , size(COEFF,1), []).';
 end%freq.
 toc

[Jan]

for j = 1 : K
   COEFF_S = COEFF_S*COEFF_S ; 
end

is not the same as

COEFF_S=COEFF_S^K;

but as

COEFF_S=COEFF_S ^ (2^K);

The matrix power is expensive and the cumulative multiplication might be faster. See FEX: mpower2

[Matt J]

How sparse is PRODOTTO3? The command

double(351,351,351)>0.7

gives an error so I don't think it could be your true code. If you really meant

    rand(351,351,351)>0.7

this will not be very sparse.

Other than that, you should avoid overhead in the ndSparse class. Do the reshaping of PRODOTTO3 prior to the loop, and convert it to a regular 2D matrix

     PRODOTTO3 = ndSparse(rand(351,351,351)>0.7);
      PD3=sparse(PRODOTTO3(:,:));

     COEFF = zeros(351,1) ; 
     COEFF(1:3) = 1 ; 
     COEFF=sparse(COEFF) ; 
     COEFF_S = sparse(size(COEFF,1),size(COEFF,1)) ; 
     tic ; 
     for l = 1:200 %freq.
     COEFF_S=reshape( COEFF.'*(PD3) , size(COEFF,1), []).';
     end%freq.
     toc

[Abi Waqas]

Thanks a lot. The real issue was the PRODOTTO3 was not very sparse. I made it more sparse then it works.

[Matt J]

After replacing with the suggested changes the main hurdle is this piece of code shown below

You could use MTIMESX( Download ). It will let you replace the loop

 COEFF_S=mtimesx(COEFF_S,A);

Answer 2

[Jan]

Why do you run the firt loop 200 times, although the loop body does not depend on the counter l? Simply omit the for loop to speed up the code by a factor 200. If the loop is required, the simplification for the forum was to heavy (or too light).

Instead of sparse matrix calculations, use logical indexing:

iter = 10;  % This is 200 in your example

S = zeros(n, n);
tic ;
n = numel(COEFF);
for l = 1:iter %freq.
   for i  = 1 : n      % Original code
      for j = 1:n
         S(i,j) = sum(COEFF .* PRODOTTO3(:,j,i));
      end
   end
end
toc

S2 = zeros(n, n);
tic ;
n = numel(COEFF);
for l = 1:iter %freq.
   for i  = 1 : n*n   % Omit the inner loop
      S2(i) = sum(COEFF .* PRODOTTO3(:,i));
   end
end
S2 = S2.';
toc

S3 = zeros(n, n);
tic ;
n = numel(COEFF);
M = (COEFF == 1);
for l = 1:iter %freq.
   for i  = 1 : n*n    % Logical indexing in a loop
      S3(i) = sum(PRODOTTO3(M,i));
   end
end
S3 = S3.';
toc

S4 = zeros(n, n);
tic;
for l = 1:iter        % Matrix multiplication
  S4 = reshape( COEFF.'*(PRODOTTO3(:,:)) , size(COEFF,1), []).';
end
toc

tic;
for l = 1:iter %freq.   % Logical indexing and matrix operation
  S5 = squeeze(sum(PRODOTTO3(COEFF == 1, :, :), 1)).';
end
toc 

 Elapsed time is 6.041851 seconds.
 Elapsed time is 5.717790 seconds.
 Elapsed time is 4.114472 seconds.
 Elapsed time is 4.665910 seconds.
 Elapsed time is 0.159271 seconds.  % !!! Logical indexing

If COEFF has more non zero entries, move the conversion to a LOGICAL out of the loop:

CC = (COEFF == 1);
for l = 1:iter %freq.
  S5 = reshape(sum(PRODOTTO3(CC, :, :), 1), n, n).';
end

 Elapsed time is 0.14 seconds.

[Abi Waqas]

Dear Jan and Matt

Thanks a lot for the brief response. The values in coefficient are not always 1 (are changing everytime big code), therefore I would like to go with S4 (matrix multiplication).

The iteration loop is there because matrices are frequency dependent which I have not shown here but 200 is the number of frequency points.

After replacing with the suggested changes the main hurdle is this piece of code shown below

if true
  for i= 1 : 250
  for j = 1 : K
  COEFF_S(:,:,i) = COEFF_S(:,:,i)*A(:,:,i) ; 
  end
end
end

K =1, it will be A(:,:,i), then K=2, it would be B(:,:,i), so forth and so on.

Meaning COEFF_S and A,B,C matrices have dimension MxN define on each frequency point which is third dimension i=200

So COEFF_S(M,N,200) and A1(M,N,200) A2(M,N,200) and AK(M,N,200)

if true

 for i= 1 : 250
    COEFF_S(:,:,i) = COEFF_S(:,:,i)*A(:,:,i) ; 
 end

 for i= 1 : 250
    COEFF_S(:,:,i) = COEFF_S(:,:,i)*B(:,:,i) ; 
 end

 for i= 1 : 250
    COEFF_S(:,:,i) = COEFF_S(:,:,i)*C(:,:,i) ; 
 end

For this, To squeeze this I have written the code like this, A is going to change and I know how to do that. But I am sure you experts have a more clever way to do that.

if true

for i= 1 : 250
    for j = 1 : K
    COEFF_S(:,:,i) = COEFF_S(:,:,i)*A_j(:,:,i) ; 
    end
end

Thank a lot for the response and help.

[Jan]

Do you mean:

 for i= 1 : 250
    COEFF_S(:,:,i) = COEFF_S(:,:,i) * A(:,:,i)  * B(:,:,i) * C(:,:,i); 
 end

Note: The definition of COEFF in the original question with only 1s and 0s was misleading. This is the problem of posting simplified code. Sorry for wasting time with this suggestion.

[Abi Waqas]

Dear Jan Simon

Thank you so much for the valuable suggestions.

It's better if we can do like this

    if true
 for i= 1 : 250
    COEFF_S(:,:,i) = COEFF_S(:,:,i)*A(:,:,i) ; 
 end

For the above code is there any faster way to do the same?

Leaving everything behind, is there no way to multiply two 3D matrices without using the for loop?

For A(M,N,i)*B(M,N,i), the only way is to do the matrix multiplication of A and B is to do for loop for i? Considering that A and B should be matrix multiplication for each value of i.

Thanks a lot for the time.

[Abi Waqas]

Dear Jan and Matt

I am still waiting for your valuable response, please.

[Matt J]

I gave you a suggestion for this here.