Matlab function for cumulative power

Is there a function in MATLAB that generates the following matrix for a given scalar r, where each row behaves somewhat like a power analog of the CUMSUM function?:

1 r r^2 r^3 ... r^n 
0 1 r   r^2 ... r^(n-1)
0 0 1   r   ... r^(n-2)
...
0 0 0   0   ... 1

[Rik]

I doubt there is a direct function. Have you tried writing one yourself?

Answer 1

[Birdman]

Try the following code. By changing r and n values, you can see the corresponding results.

r=9;n=4;
A=zeros(n+1,n+1);
for i=1:size(A,1)
    for j=1:size(A,2)
        if (j-i)<0
            A(i,j)=0;
        else
            A(i,j)=r^(j-i);
        end
    end
end
A

[Herr K]

Thanks. This is straightforward.

Answer 2

[Rik]

This code does what you ask without loops.

%define inputs
r=9;
n=4;

[a,b]=meshgrid(0:n);
exponents=a-b;
exponents(exponents<0)=NaN;
result=r.^exponents;
result(isnan(result))=0;

[Birdman]

This code is slower than mine although you avoided nested for loops.

[Rik]

Strange. I did a tic,toc to check and mine was about twice as fast, but with the code below yours is 10 times faster. It probably also depends on n and r.

r=9;n=4;
timeit(@() option_loop(r,n))
timeit(@() option_grid(r,n))

function A=option_loop(r,n)
A=zeros(n+1,n+1);
for i=1:size(A,1)
    for j=1:size(A,2)
        if (j-i)<0
            A(i,j)=0;
        else
            A(i,j)=r^(j-i);
        end
    end
end
end
function result=option_grid(r,n)
[a,b]=meshgrid(0:n);
exponents=a-b;
exponents(exponents<0)=NaN;
result=r.^exponents;
result(isnan(result))=0;
end

[Birdman]

Your answer is neat but my point is there is no reason to avoid nested for loops because there is no dramatic time consumption.

[Rik]

This is actually a nice illustration of the fact that a non-loop version isn't always faster. In this case (at least on my computer with Windows 10 and R2019a) the looped version is faster up to about n=30. For huge values of n there may very well be a tangible benefit (or if this code is going to be run very often).

clc,clear
r=9;
n_list=[1:100 200:100:1000];
t=zeros(2,numel(n_list));
for it=1:size(t,2)
    n=n_list(it);
    t(1,it)=timeit(@() option_loop(r,n));
    t(2,it)=timeit(@() option_grid(r,n));
end

figure(1),clf(1)
plot(n_list,t(1,:),n_list,t(2,:))
legend({'loop','grid'})
xlabel('n'),ylabel('time')

[Herr K]

Thanks for your answer. Now I learn the meshgrid function as well!