Using vectorization on a problem to reduce For loops, how to use conditional statement?

I have the following function code and I want to reduce the number of for loops it has:

function [y, zf] = MyFunction(x, b, zi)

y = zeros([length(x) 1]);

for n = 1:length(x)
  for k=1:length(zi)
       if n<k
           y(n) = y(n) + b(k)*zi(length(zi)+(n-k)+1);
       else
           y(n) = y(n) + b(k)*x(n-k+1);
       end
   end
end

zf = x(length(x)-length(zi)+1:length(x));

I manage to do kindo of do that with vectorization, but I am having problems with the if statement, how could I solve it?

function [y, zf] = MyFunction(x, b, zi)

y = zeros([length(x) 1]);

n=1:1:length(x); % use of vectorization
for k=1:length(zi)
   if n<k % problem with if
       y = y + b(k)*zi(length(zi)+(n-k)+1);
   else
       y = y + b(k)*x(n-k+1);
   end
end

zf = x(length(x)-length(zi)+1:length(x));

[Rik]

What you could try is one step further: use logical indexing. That way you can separate the two cases.

[Pedro Vieira]

I read this documentation about Indexing. So if I do something like:

if n(1:length(x))<k

It should work?

Answer 1

[Rik]

You can find an explanation of logical indexing on that page you linked.

The most important thing if you want to vectorize an operation is that you need to remove dependency on other operations. I have done that in the code below. It is quite clunky and could be further optimized, but it does not require loops. Timing it with small example inputs shows how good Matlab is in optimizing loops (run a few times and you'll see the difference).

I don't fully understand what your code is supposed to do, so I don't really see an easy way to optimize this code.

clearvars,clc%only use clearvars for debugging
x=rand(10,1);
b=rand(6,1);zi=rand(size(b));

tic
[y, zf] = MyFunction(x, b, zi);
toc

tic
%alternative:
[n,k]=ndgrid(1:length(x),1:length(zi));
y1=zeros(size(n));
L=n<k;
y1(L) = y1(L) + b(k(L)).*zi(length(zi)+(n(L)-k(L))+1);
L=~L;%flip to do the else
y1(L)=y1(L) + b(k(L)).*x(n(L)-k(L)+1);
y1=sum(y1,2);%sum over second dimension to re-apply the recursive behavior
toc

fprintf('this number should be (almost) 0: %.5f\n',max(abs(y-y1)))

function [y, zf] = MyFunction(x, b, zi)
y = zeros([length(x) 1]);

for n = 1:length(x)
  for k=1:length(zi)
       if n<k
           y(n) = y(n) + b(k)*zi(length(zi)+(n-k)+1);
       else
           y(n) = y(n) + b(k)*x(n-k+1);
       end
   end
end
zf = x(length(x)-length(zi)+1:length(x));
end

[Rik]

Did this suggestion solve your problem? If so, please consider marking it as accepted answer. It will make it easier for other people with the same question to find an answer. If this didn't solve your question, please comment with what problems you are still having.