Compute double sum in matlab efficiently?

Question

I am looking for an optimal way to program this summation ratio. As input I have two vectors v_mn and x_mn with (M*N)x1 elements each.

The ratio is of the form:

The vector x_mn is 0-1 vector so when x_mn=1, the ration is r given above and when x_mn=0 the ratio is 0.

The vector v_mn is a vector which contain real numbers.

I did the denominator like this but it takes a lot of times.

function r_ij = denominator(v_mn, M, N, i, j)
%here x_ij=1, to get r_ij.
S = [];
for m = 1:M
  for n = 1:N
    if (m ~= i)
      if (n ~= j)
        S = [S v_mn(i, n)];
      else
        S = [S 0];
      end
    else
      S = [S 0];
    end
  end
end
r_ij = 1+S;
end

Can you give a good way to do it in matlab. You can ignore the ratio and give me the denominator which is more complicated.

EDIT: I am sorry I did not write it very good. The i and j are some numbers between 1..M and 1..N respectively. As you can see, the ratio r is many values (M*N values). So I calculated only the value i and j. More precisely, I supposed x_ij=1. Also, I convert the vectors v_mn into a matrix that's why I use double index.

Answer 1

If you reshape your data, your summation is just a repeated matrix/vector multiplication.

Here's an implementation for a single m and n, along with a simple speed/equality test:

clc

%# some arbitrary test parameters
M = 250;
N = 1000;
v = rand(M,N);   %# (you call it v_mn)
x = rand(M,N);   %# (you call it x_mn)

m0 = randi(M,1); %# m of interest
n0 = randi(N,1); %# n of interest 


%# "Naive" version
tic
S1 = 0;
for mm = 1:M %# (you call this m')
    if mm == m0, continue; end
    for nn = 1:N %# (you call this n')
        if nn == n0, continue; end
        S1 = S1 + v(m0,nn) * x(mm,nn);
    end
end
r1 = v(m0,n0)*x(m0,n0) / (1+S1);
toc


%# MATLAB version: use matrix multiplication!
tic

ninds = [1:m0-1 m0+1:M];
minds = [1:n0-1 n0+1:N];
S2 = sum( x(minds, ninds) * v(m0, ninds).' );
r2 = v(m0,n0)*x(m0,n0) / (1+S2);

toc


%# Test if values are equal
abs(r1-r2) < 1e-12

Outputs on my machine:

Elapsed time is 0.327004 seconds.   %# loop-version
Elapsed time is 0.002455 seconds.   %# version with matrix multiplication
ans =  
     1                              %# and yes, both are equal

So the speedup is ~133×

Now that's for a single value of m and n. To do this for all values of m and n, you can use an (optimized) double loop around it:

r = zeros(M,N);
for m0 = 1:M   
    xx = x([1:m0-1 m0+1:M], :);
    vv = v(m0,:).';
    for n0 = 1:N
        ninds    = [1:n0-1 n0+1:N];        
        denom    = 1 + sum( xx(:,ninds) * vv(ninds) );
        r(m0,n0) = v(m0,n0)*x(m0,n0)/denom;        
    end
end

which completes in ~15 seconds on my PC for M = 250, N= 1000 (R2010a).

EDIT: actually, with a little more thought, I was able to reduce it all down to this:

denom = zeros(M,N);
for mm = 1:M    
    xx = x([1:mm-1 mm+1:M],:);
    denom(mm,:) = sum( xx*v(mm,:).' ) - sum( bsxfun(@times, xx, v(mm,:)) );    
end
denom = denom + 1;

r_mn = x.*v./denom;

which completes in less than 1 second for N = 250 and M = 1000 :)