How to sum over indices i+j=k without using a for loop?

I have an matrix β and a vector n. I am trying to compute two sums:

1. 
2. 

(k is an integer). So I should have answers for and . Is there a way to compute these sums without using for loops or computationally expensive Matlab matrix manipulation functions? The first sum in particular is giving me a headache. Thanks!

[Jan]

Why do you want to avoid a loop? This wish occurs frequently in questions in this forum, but I do not see why? Since Matlab 6.5 the JIT accelerates loops efficiently and the old rumor "do not use loops in Matlab" is outdated for 20 years now.

Did you implement this with loops already? Then please post the corresponding code.

What are typical values for m?

[Jackie Taylor]

These are sums in the middle of a simulation, so I'm already looping through time and space. and β needs to be re-calculated at every grid cell, for every timestep. I tried to tease out only the relevant bits of code from a much larger script, copied below.

% spatial grid
delz = .1; % m
zmax = 6; % m
z = 0.5*delz:delz:zmax-0.5*delz; % spatial grid
p = length(z);

% time steps
delt = 20; % s
numDays = 7; % how many days to model
t = 0:delt:60*60*24*numDays; % simulation time in seconds
m = length(t);

% volume
delv = 1; % cell
kmax = 501; 
vol = 1:delv:vmax; 
r = length(vol); 

% initialize sums
qin = NaN(p,r);
qout = NaN(p,r);

for j = 1:m-1 % time loop

% initialize sums
qin = zeros(p,r);
qout = zeros(p,r);

    for i = 1:p % space loop           
        beta = ...; % rxr matrix
        for k = 1:kmax
            for s = 1:k-1          
                qin(i,k) = qin(i,k) + beta(s,k-s)*nold(i,s)*nold(i,k-s); % nold is pxr matrix, redefined every timestep
            end

            for s = 1:r-k                
                qout(i,k) = qout(i,k) + beta(s,k)*nold(i,s)*nold(i,k);
            end
    end    
end

[Jan]

@Jackie Taylor: What is vmax and nold? I try:

vol = 1:delv:kmax;  % Instead of vmax
n   = rand(p, r);   % nold -> n

[Jackie Taylor]

Good catch, vmax = kmax. I'm going to try your suggested inner loop edits and get back to you on the results. Thanks for all the help!

Answer 1

[Jan]

Because beta is determined inside the for i and for j loops, we can improve only the inner loops:

r    = 501;
kmax = 501;
p    = 60;
beta = rand(r, r);
n    = rand(p, r);
qin  = zeros(p, r);
qout = zeros(p, r);
i    = 1;
tic
for rep = 1:100  % Only to get a more accurate time measurement, remove in real code!
    for k = 1:kmax
        for s = 1:k-1          
            qin(i,k) = qin(i,k) + beta(s,k-s) * n(i,s) * n(i,k-s);
        end
        for s = 1:r-k                
            qout(i,k) = qout(i,k) + beta(s,k) * n(i,s) * n(i,k);
        end
    end
end
toc

Elapsed time is 0.449715 seconds.

qin  = zeros(p, r);
qout = zeros(p, r);
tic
ni = n(i, :);
for rep = 1:100
    for k = 1:kmax
        a = 0;
        for s = 1:k-1          
            a = a + beta(s, k-s) * ni(s) * ni(k-s);
        end
        qin(i, k) = qin(i, k) + a;
        
        a = 0;
        for s = 1:r-k                
            a = a + beta(s, k) * ni(s);
        end
        qout(i, k) = qout(i, k) + a * ni(k);
    end
end
toc

Elapsed time is 0.220807 seconds.

Almost 2 times faster with just avoiding some indexing. Because this is a part of the problem only the total speedup will be smaler. My trials to vectorize this are about 10 times slower. In addition the timings measure in the online version differ from a local version. So please check and post, if the modified version is faster at all.

[Jackie Taylor]

This works much better than what I had in mind. Thanks!

Answer 2

[Walter Roberson]

If you have i+j = k, then you start at row k-1, column 1 as the first point (if k >= the number of rows). The next point would be row k-2, column 2. The third point would be row k-3, column 3. And so on.

In terms of linear indices, the first one is (k-1) + rows*(1-1); the second one is (k-2)+rows*(2-1), the third is (k-3)+rows*(3-1) and so on. Those are k-1, k-2+rows, k-3+2*rows, k-4+3*rows and so on. The difference between those is [rows-1], [rows-1], [rows-1], [rows-1] and so on.

Therefore the linear indices are k-1:rows-1:SOME_ENDPOINT . And that can be implemented as a vector without using a loop.

You will need some logic for the case where k is less than the number of rows or columns.