Iteration without any loop.

I am working on B-spline, this code is to get a matrix as an output.

function [M] = Iterations(N) % N is an input variable
I=eye(N+1);M=zeros(N+1,N+1); % I is a diagnal matrix and M is inicially a zero matrix
for i=0:N/2 % interations 
n=i:N-i;
b=(-1).^(N-n-i).*factorial(N-i)./factorial(n)./factorial(N-n-i).*(factorial(N)./factorial(N-i)./factorial(i));
M(i+1,i+1:N+1-i)=b;
end
M = rot90(fliplr(M),1)+M-diag(M).*I; % Output matrix
end
Ans 
M =

    -1     3    -3     1
     3    -6     3     0
    -3     3     0     0
     1     0     0     0

I was wondering is there anyway i can remove a loop from inside the function here?

Answer 1

[Walter Roberson]

[M] = Iterations(3)

M = 4×4

-1 3 -3 1 3 -6 3 0 -3 3 0 0 1 0 0 0

function [M] = Iterations(N) % N is an input variable
I=eye(N+1);M=zeros(N+1,N+1); % I is a diagnal matrix and M is inicially a zero matrix
i=0:N/2; % interations 
n=(0:N).';
b = (-1).^(N-n-i).*gamma(N-i+1)./factorial(n)./gamma(N-n-i+1).*(factorial(N)./gamma(N-i+1)./factorial(i));
M = [tril(b).'; zeros(size(b.'))];
M = rot90(fliplr(M),1)+M-diag(M).*I;
end

[Jan]

This fails for even N:

Matrix dimensions must agree.
Error in kjh>Iterations_4 (line 84)
  M = rot90(fliplr(M),1)+M-diag(M).*I;

Answer 2

[Jan]

Why do you want to do this without a loop? When the purpose is an acceleration, omit the expensive factorial and power operations at first:

function M = Iterations(N)
F  = [1, cumprod(1:21)];   % Emulated factorial()
sF = F;
sF(2:2:22) = -sF(2:2:22);  % Include expensive power operation -1^(N-n-i)
FN = F(N+1);
M  = zeros(N+1, N+1);
for i = 0:N/2
   t = N - i + 1;
   b = F(t) ./ F(i+1:t) ./ sF(t-i:-1:1) .* (FN ./ F(t) ./ F(i+1));
   M(i+1, i+1:t) = b;
   M(i+1:t, i+1) = b;
end
end

For N = 10 this is needs 16% of the runtime of the original.

But due to the factorial the output will be accurate for N <= 21. If you need such a small set of inputs only, it would be more efficient, to create a list of possible values once and store it in a persistent varibale.

function M = IterationsX(N)
persistent List
if isempty(List)
    List = cell(1, 21);
    for k = 1:21
        List{k} = Iterations(k);
    end
end
M = List{N};
end

[hyder abedi]

Thanks for improving the code and giving me a different way of writing it. I was just curious whether this or any iteration can be done without a loop at all. yet the question remains.

And how did you calculate the runtime?

[Image Analyst]

You can use tic and toc to measure elasped seconds.

For loops are not always slow. How big is your N? Unless your N is like hundreds of millions, I doubt the loop will take more than a fraction of a second. For N = 10 or 21, you're looking at an elapsed time of like microseconds.

[hyder abedi]

It depends on number of variables in my dataset, right now its 8 but it can be 100 or more.

[Jan]

factorial(100) is 9.3e157. BSplines with so many points cannot be calculated accurately anymore, because the double values contain the first 16 digits only.

[Walter Roberson]

When you have factorial divided by factorial, consider using pochhammer; https://www.mathworks.com/help/symbolic/pochhammer.html

[hyder abedi]

Thanks again Jan i wasn't aware that double values contain the first 16 digits only.

[hyder abedi]

Thanks Walter Roberson.