Thread Subject: is there any way to vectorize the following code with two loops?

Dear all,

It is the second time for me to post this thread. I still haven't get the solution yet. I hope the experts here can answer my question. As I heard, most loops can be vectorized to speed up the program. Is there any way to accelerate my current code? If this is already the efficient way using Matlab to solve this question, I have to turn to other languages, say, Fortran, to speed it. That means I need to spend more time to learn a new language for me. So, please, help me!

The original post:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
I have written a code as below, using two loops. Is there a smart way to use vectorization instead of using loops to speed it? I mean, without loops or only one loop? Since the current code is inside some main program and it does substantially impose the computation burden when I run the whole program each time.

The code is to solve the following problem:
to express a matrix D_{ni} where the subscript {ni} means the nth row, the ith column element in matrix D. D is a N*N matrix. D must satisfy the following condition (the writing follows the format of Latex):
D_{ni}=\sum_{j=1}^{N} y_{nij} * b_{j} * {\min_{m=1,2, ..., N} [c_{m}*K_{nm}*L_{mj}]}

where D, K, L are N*N matrix. b and c are N*1 vector. y is N*N*N (3-dimension). Suppose b, c, K, L are given. We know y_{nij}=0 if i~=m and y_{nij}=1 if i==m. I mean, if i is the index which minimizes the last part (inside the \min symbol), then this y_{nij}=1. If there are more than one index to minimize the last part, we evenly give the weight to each index. Therefore, ultimately, we have \sum_{i} y_{nij}=1
i, j, n, m are all indexes which is taken from 1 to N. Now we need to express this D matrix. We can know that we only need to consider when i (in D_{ni}) equals the index m (i.e., the index which can minimize the last part) and all other D_{ni}=0

My code is as below.
y=zeros(N,N,N);
D=zeros(N);
for n=1:N
    for j=1:N
        val=c.*K(n,:)'.*L(:,j);
        minval=min(val);
        minj=find(val==minval);
        y(n,minj,j)=1/length(minj);
        D(n,minj)=D(n,minj)+minval*b(j)/length(minj);
    end
end

If the two loops above can be vectorized, that would be a great help to my program! I hope someone can help me out. Many thanks in advance!
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

BTW, it doesn't matter whether or not y shows up in the expression of D since all other y =0 and only y_{n,minj,j) matters. When constructing D, we only need to construct D(n,minj) that's why y does not appear there. Please see below.
        y(n,minj,j)=1/length(minj);
        D(n,minj)=D(n,minj)+minval*b(j)/length(minj);
Or you can write as this: D(n,minj)=D(n,minj)+minval*b(j)*y(n,minj,j);


All my best,
Yao (Amber)

This is the vectorized code:

CK = bsxfun(@times,c',K);
CKL = bsxfun(@times,CK,reshape(L,[1 N N]));
clear CK
minval = min(CKL,[],2);
mineq = bsxfun(@eq,CKL,minval);
clear CKL
il = 1./sum(mineq,2);
y = bsxfun(@times, mineq, il);
D = bsxfun(@times, minval, il);
D = bsxfun(@times, D, reshape(b, [1 1 N]));
D = bsxfun(@times, D, mineq);
D = sum(D,3);

% Bruno

Hi Bruno,

Thank you so much! I have run the program, and it does increase the speed when N is in a reasonable size (N=30, N=130, etc.)! However, it's interesting that if N is very large (for e.g., I tried N=330), the vectorization code is not as efficient as the previous loops. But I expect vectorization will be more efficient when N is growing larger and larger. Am I right? Or something missed here?

I feel I've so much to learn. Thank you a lot!!

Cheers,
Yao (Amber)

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h770om$d6s$1@fred.mathworks.com>...
> This is the vectorized code:
>
> CK = bsxfun(@times,c',K);
> CKL = bsxfun(@times,CK,reshape(L,[1 N N]));
> clear CK
> minval = min(CKL,[],2);
> mineq = bsxfun(@eq,CKL,minval);
> clear CKL
> il = 1./sum(mineq,2);
> y = bsxfun(@times, mineq, il);
> D = bsxfun(@times, minval, il);
> D = bsxfun(@times, D, reshape(b, [1 1 N]));
> D = bsxfun(@times, D, mineq);
> D = sum(D,3);
>
> % Bruno

The benefits of vectorization do depend on the size of the operands. For example, for N=300, this is faster than Bruno's vectorized code AND your loops on my machine:



y = zeros(N,N,N,'double');
D = zeros(N,'double');
b = b.';
for n=1:N
    V = c.*K(n,:)';
    val = bsxfun(@times,V,L);
    minval = min(val);
    M = minval.*b;
    for jj=1:N
        minj = (val(:,jj)==minval(jj));
        NNZ = sum(minj);
        y(n,minj,jj) = 1/NNZ;
        D(n,minj) = D(n,minj) + M(jj)/NNZ;
    end
end

Thank you very much! I will try this too :)

"Matt Fig" <spamanon@yahoo.com> wrote in message <h7733n$cbm$1@fred.mathworks.com>...
> The benefits of vectorization do depend on the size of the operands. For example, for N=300, this is faster than Bruno's vectorized code AND your loops on my machine:
>
>
>
> y = zeros(N,N,N,'double');
> D = zeros(N,'double');
> b = b.';
> for n=1:N
> V = c.*K(n,:)';
> val = bsxfun(@times,V,L);
> minval = min(val);
> M = minval.*b;
> for jj=1:N
> minj = (val(:,jj)==minval(jj));
> NNZ = sum(minj);
> y(n,minj,jj) = 1/NNZ;
> D(n,minj) = D(n,minj) + M(jj)/NNZ;
> end
> end

"Yao Li" <yaoamber@gmail.com> wrote in message <h772od$jip$1@fred.mathworks.com>...
> Hi Bruno,
>
> Thank you so much! I have run the program, and it does increase the speed when N is in a reasonable size (N=30, N=130, etc.)! However, it's interesting that if N is very large (for e.g., I tried N=330), the vectorization code is not as efficient as the previous loops. But I expect vectorization will be more efficient when N is growing larger and larger. Am I right? Or something missed here?
>

The vectorized code is slow because the mininum draw-equality is treated by expansion the third dimensional array, which is a costly strategy.

Bruno