MATLAB: Getting rid of nesting for loops which rely on reordering

matrixmatrix manipulationnested for loop

I have two K x N matrices B and H.. I need to calculate a new matrix z according to:
     z(k,n) = a + B(k,n)*H(k,n) + ...
             sum{j = 1...(i-1)} 2B(j,n)*H(j,n)t +...
             sum{j = (i+1)...K)} B(j,n)*H(j,n)
where a and t are scalar constants and i depends on the "rank" in the column if H is sorted descending, i.e. for
     H = [1,2;5,6;3,4];
Column ranks are 2,3,1 for both columns. Sorry if this is confusing, I'm doing my best.
Anyways, I need to do this calculation for every element of z and I'm always only worried about the rank in each column and all terms in z are of elements of H and B in the same column.
My naive attempt to solve is using nested for loops:
    z = H.*B + a;
    for n=1:N       
        [Hsorted, idx] = sort(H(:,n),'descend');  % look at each column      
        for i=1:K                   %for each element of z(:,n)
            for j=1:K               %look at every other element in H(:,n)
                if j ~= i 
                    if  j < i                             
                        z(idx(i),n) = z(idx(i),n) + 2*B(idx(j),n)*Hsorted(j)*t;
                    else                               
                        z(idx(i),n) = z(idx(i),n) + B(idx(j),n)*Hsorted(j);
                    end
                end
            end  
        end 
    end
As you can probably imagine, this takes a lot of time for large K and N. I also am doing this as part of an iterative algorithm and a calculation like this gets done 2 or 3 times on each iteration. So, even for pretty small K and N, the time really adds up.
I've figured out that I can sort H, then reorder B according to the order of H with
     [Hsorted,orderH] = sort(H,'descend');
     [m,n] = size(B);
     Bsorted = B(sub2ind([m,n],orderH,repmat(1:n,m,1)));
and I'm guessing (hoping) that I can somehow replicate the ranked summation part of calculating the elements of z with matrix operations. This is where I'm stuck. Any help/advice or even a direction would be really appreciated

Best Answer

      [K,N]=size(H);
       z = H.*B + a;
       [Hsorted, orderH] = sort(H,1,'descend');
       %sort B consistently with Hsorted
       idx=sub2ind([K,N],orderH,repmat(1:N,K,1));
       Bsorted=B(idx);
           e=ones(K);
           W=(triu(e,1)+2*t*tril(e,-1));
       z(idx)=z(idx)+W*(Bsorted.*Hsorted);
Related Question