Matlab double sum over vectors

I want to calculate , in which k and k' are vectors with x and y components. Both k and k' lie in the 1st Brillouin zone: -pi<kx<pi, -pi<ky<pi, -pi<kx'<pi, -pi<ky'<pi. How to write Matlab code to perform the sum over k and k'? My code doesn't work:

sum = 0;
kvec = linspace(-pi,pi,N); 

for j1=1:1:length(kvec)
    kx = kvec(j1);
    for j2=1:1:length(kvec)
        ky = kvec(j2);
        for j3=1:1:length(kvec)
            kxprime = kvec(j3);
            for j4=1:1:length(kvec)
                kyprime = kvec(j4);
                
                sum = sum + f(kx,ky,kxprime,kyprime)
            end
        end
    end
end

[Walter Roberson]

it depends on whether f is vectorized

[CC SS]

f is a scalar in terms of kx, ky, kx' and ky', not a vector.

[Dyuman Joshi]

Can you share the definition of f?

[CC SS]

f = 2* (cos(kx) + cos(ky)) - 0.5* (cos(kxprime) + cos(kyprime)).

Answer 1

[Walter Roberson]

Different versions, and their timings.

Except... in different runs, the timings especially for the first version varied by more than a factor of 10, so the values shown here for a single run might not be representative of real use.

format long g

tic
N = 10;
kvec = linspace(-pi,pi,N);
ck = cos(kvec);
ca = reshape(ck, [], 1);
cb = reshape(ck, 1, []);
cc = reshape(ck, 1, 1, []);
cd = reshape(ck, 1, 1, 1, []);
f = (ca + cb)*2 - (cc + cd)/2;
out1 = sum(f,'all')

out1 = 
         -2999.99999999997

toc

Elapsed time is 0.022545 seconds.

tic
f = @(a,b,c,d) 2*(cos(a) + cos(b)) - 0.5*(cos(c) + cos(d));
%Random value for N for example
N = 10;
kvec = linspace(-pi,pi,N); 
[kx,ky,kxprime,kyprime] = ndgrid(kvec);
arr = f(kx,ky,kxprime,kyprime);
out2 = sum(arr,'all')

out2 = 
         -2999.99999999997

toc

Elapsed time is 0.007611 seconds.

tic
out3 = 0;
N = 10;
kvec = linspace(-pi,pi,N);
f = @(k,kprime) 2* (cos(k(1)) + cos(k(2))) - 0.5* (cos(kprime(1)) + cos(kprime(2)));
for j1=1:1:length(kvec)
    kx = kvec(j1);
    for j2=1:1:length(kvec)
        ky = kvec(j2);
        for j3=1:1:length(kvec)
            kxprime = kvec(j3);
            for j4=1:1:length(kvec)
                kyprime = kvec(j4);                
                out3 = out3 + f([kx,ky],[kxprime,kyprime]);
            end
        end
    end
end
out3

out3 = 
         -3000.00000000004

toc

Elapsed time is 0.011396 seconds.

tic
N = 10;
kvec = linspace(-pi,pi,N);
ck = cos(kvec);
ca = reshape(ck, [], 1);
cb = reshape(ck, 1, []);
cc = reshape(ck, 1, 1, []);
cd = reshape(ck, 1, 1, 1, []);
f = bsxfun(@plus, bsxfun(@plus, ca, cb), bsxfun(@plus, cc, cd));
out4 = sum(f(:))

out4 = 
         -3999.99999999997

toc

Elapsed time is 0.010282 seconds.

Answer 2

[Torsten]

Note that -pi<=kx<=pi, -pi<=ky<=pi, -pi<=kx'<=pi, -pi<=ky'<=pi because of your linspace choice.

And most probably you need to normalize the sum somehow because at the moment, it depends strongly on N.

sum = 0;
N = 10;
kvec = linspace(-pi,pi,N);
f = @(k,kprime) 2* (cos(k(1)) + cos(k(2))) - 0.5* (cos(kprime(1)) + cos(kprime(2)));
for j1=1:1:length(kvec)
    kx = kvec(j1);
    for j2=1:1:length(kvec)
        ky = kvec(j2);
        for j3=1:1:length(kvec)
            kxprime = kvec(j3);
            for j4=1:1:length(kvec)
                kyprime = kvec(j4);                
                sum = sum + f([kx,ky],[kxprime,kyprime]);
            end
        end
    end
end
sum

sum = -3.0000e+03

Answer 3

[Dyuman Joshi]

Using 4 nested for loops will be take quite a good amount of time to run, specially if N is a comparetively big value.

You can vectorize your code -

f = @(a,b,c,d) 2*(cos(a) + cos(b)) - 0.5*(cos(c) + cos(d));

%Random value for N for example
N = 10;

kvec = linspace(-pi,pi,N); 
[kx,ky,kxprime,kyprime] = ndgrid(kvec);

arr = f(kx,ky,kxprime,kyprime);

Also, it's not a good idea to use inbuilt function names as variables, in your case - sum

out = sum(arr,'all')

out = -3.0000e+03