Efficient way to multiply a vector of size (Nx,1) with a matrix of size (Nx+1,Nx+1)

My code repeatedly computes a vector of size (Nx,1) by a matrix of size (Nx+1,Nx+1):

clc; clear all;
Nx = 32;
a = 1.5;
b = 2;
gamma = rand(Nx,1); 
D = rand(Nx+1,Nx+1);
D2 = D*D;
identy = eye(Nx+1,Nx+1);

% Computation performed by my code millions of times
for j=1:Nx
  A = a*D2 + b*D + gamma(j)*identy;
  % Can gamma(j)*identy be avoided through vectorization or a faster
  % technique?
end

Is there an alternative way to compute A without using the for loop (would vectorization or repmat also work)?

[Matt J]

In your example code, A does not evolve throughout the loop. Rather it is repeatedly over-written resulting in

 A = a*D2 + b*D + gamma(Nx);

This is probably not what you want, but the reader cannot tell what it should really be doing.

[Matthew Kehoe]

I agree @Matt J. I didn't want to ask about all of my code (as it might have appeared to be a difficult question) so I only asked about a small portion of my code. This might have been a mistake. The answer provided by the cyclist has helped me fix/improve my own code significantly.

Answer 1

[the cyclist]

If you permute gamma to be a vector along the dimension 3, then you can multiply it by identy, and the result will be a 33x33x32 array, where each "slice" along dimension 3 is the multiple with the corresponding value of gamma.

So, this calculation will do all of the calculations of your for loop (and store them all, rather than overwriting as your code does).

A = a*D2 + b*D + permute(gamma,[3,2,1]) .* identy;

You may now have a memory problem, though, if your arrays are large.

[Matthew Kehoe]

I see, I think that it depends on how you manipulate A inside the for loop. If you have other actions inside the for loop then computing A inside the for loop appears to be faster.

clc; clear all;
Nx = 32;
a = 1.5;
b = 2;
gamma = rand(Nx,1); 
D = rand(Nx+1,Nx+1);
D2 = D*D;
D_start = D(1,:);
D_end = D(end,:);
identy = eye(Nx+1,Nx+1);
n_min = 1.2;
d_min = 1.6;
n_max = 7;
d_max = 1.0;


tic
% Computation performed by my code millions of times
% Run 1000 computations
for ntests=1:1000
  for j=1:Nx
    A = a*D2 + b*D + gamma(j)*identy;
    A(end,:) = n_min*D_end;
    A(end) = A(end) + d_min;

    A(1,:) = n_max*D_start;
    A(1,1) = A(1,1) + d_max;
  end
end
toc
% Elapsed time is 0.042577 seconds.


tic
% Move A outside of the j loop
A2 = a*D2 + b*D + permute(gamma,[3,2,1]) .* identy;
% Computation performed by my code millions of times
% Run 1000 computations
for ntests=1:1000
  for j=1:Nx
    A2(end,:,j) = n_min*D_end;
    A2(end,end,j) = A2(end,end,j) + d_min;

    A2(1,:,j) = n_max*D_start;
    A2(1,1,j) = A2(1,1,j) + d_max;
    A_orig = A2(:,:,j);
  end
end
toc
% Elapsed time is 0.124658 seconds.

Therefore, if you need to change indices then it might be necessary to do everything inside the for loop.

[the cyclist]

But those other operations don't need to be done in the loop, either. When I ran this on my machine, the fully vectorized code is about 3 times faster.

NTEST = 10000;

Nx = 32;
a = 1.5;
b = 2;
gamma = rand(Nx,1); 
D = rand(Nx+1,Nx+1);
D2 = D*D;
D_start = D(1,:);
D_end = D(end,:);
identy = eye(Nx+1,Nx+1);
n_min = 1.2;
d_min = 1.6;
n_max = 7;
d_max = 1.0;

tic
% Computation performed by my code millions of times
% Run 1000 computations
for ntests=1:NTEST
  for j=1:Nx
    A = a*D2 + b*D + gamma(j)*identy;
    A(end,:) = n_min*D_end;
    A(end) = A(end) + d_min;

    A(1,:) = n_max*D_start;
    A(1,1) = A(1,1) + d_max;
  end
end
toc

Elapsed time is 0.311290 seconds.

tic
% Move A outside of the j loop
A2 = a*D2 + b*D + permute(gamma,[3,2,1]) .* identy;
% Computation performed by my code millions of times
% Run 1000 computations
for ntests=1:NTEST
  for j=1:Nx
    A2(end,:,j) = n_min*D_end;
    A2(end,end,j) = A2(end,end,j) + d_min;

    A2(1,:,j) = n_max*D_start;
    A2(1,1,j) = A2(1,1,j) + d_max;
    A_orig = A2(:,:,j);
  end
end
toc

Elapsed time is 1.313793 seconds.

tic
% Move A outside of the j loop
A3 = a*D2 + b*D + permute(gamma,[3,2,1]) .* identy;
% Computation performed by my code millions of times
% Run 1000 computations
for ntests=1:NTEST
    A3(end,:,:) = repmat(n_min*D_end,[1,1,Nx]);
    A3(end,end,:) = A3(end,end,:) + d_min;
    A3(1,:,:) = repmat(n_max*D_start,[1,1,Nx]);
    A3(1,1,:) = A3(1,1,:) + d_max;
end
toc

Elapsed time is 0.108381 seconds.

[the cyclist]

I just noticed that in your A2 code, you moved the initial calculation of A2 not only outside the "j" loop, but also out of the "ntests" timing loop (presumably inadvertently). Moving that back inside the timing loop slows things down. But all these timings will likely depend on the sizes of your arrays, and there may be further optimizations.

[Matthew Kehoe]

Yes, I forgot to put A2 inside the ntests loop. Thanks for all of your help @the cyclist. Would it be alright if I asked you one more question about my code (I almost removed the j loop completely but got stuck on one part)? My code has one more place that I can optimize. Would it be alright if I asked you this in a comment (as opposed to making a new question)?

[the cyclist]

That's fine, but if the question is self-contained, it might be better to open a new one, because it might get wider exposure.

You can always tag me as you did in your comment, and I will get a notification.