Star
Star

Reputation: 2299

How to vectorize code with 3 loops that generates a matrix

Could you help me to vectorise this Matlav code constructing a matrix A of dimension MNx(2+N-1)xR in order to speed it up? At the moment it takes approx. 8 sec

INITIALISING

R=200;
M=400;
N=20;
B=[kron((1:1:M)', ones(N,1)) repmat((1:1:N)', M,1)]; %(MN)x(2)

B looks like

B=[1 1;
   1 2;
   ...;
   1 20;
   2 1;
   2 2;
   ...
   2 20;
   ...
   400 20]

CODE

A=[ repmat(B,1,1,R)  zeros(M*N,N-1,R)]; %Allocate the space
                                        %(MN)x(2+N-1)x(R)
                                        % I want A(:,1:2,r)=B for r=1,...,R


%Fill the remaining columns of A in the following way
for r=1:R
    utemp1=randn(M*N, N-1); %Generate a matrix of random 
                            %numbers of dimension (M*N)x(N-1)
    utemp2=randn(M, N);  %Generate a matrix of random 
                         %numbers of dimension (M)x(N)
    utemp3=zeros(M*N,N-1); %Generate a matrix of random 
                           %numbers of dimension(M)x(N-1)
    for m=1:M
        for j=1:N
            utemp3((m-1)*N+j,:)= utemp2(m,j)+[utemp2(m,1:j-1) utemp2(m,j+1:N)]; %(1)x(N-1)
              %e.g. if m=2, j=3, I want to fill the 23th row of utemp3
              %with utemp2(2,3)+[utemp2(2,1:2) utemp2(m,4:20)];

              %e.g. if m=4, j=1, I want to fill the 61st row of utemp3
              %with utemp2(4,1)+[utemp2(4,2:20)];
        end
    end
    A(:,3:end,r)=utemp1+utemp3;     %sum utemp1
end

SOME EXPLANATIONS

for r=1,...,R

A is such that

for m=1,...,M and for j=1,...,N

the row in A(:,:,r) starting with [m j] in the first two columns is filled in the remaining (N-1) columns with uj+uh+ujh for each h~=j, where uj, uh, ujh are i.i.d standard Gaussian numbers that can be found in utemp1 and utemp2.

Upvotes: 0

Views: 46

Answers (1)

rahnema1
rahnema1

Reputation: 15867

You can precompute the indices and use them in 500 iterations:

idx = repmat(reshape(1:M*N,M,N).',1,N);
idx = reshape(idx(logical(kron(~eye(N),ones(1,M)))),N-1,[]).';
for k=1:500
    for r=1:R
        utemp2=randn(M*N,1);
        A(:,3:end,r)=randn(M*N, N-1)+bsxfun(@plus,utemp2,utemp2(idx) );
    end
end

However allocating large matrices ,specially of repeated elements, and vectorizing operations on it is not always the most efficient way. There are built-in functions that directly operate on the original array avoiding repeating elements of the array.

Upvotes: 1

Related Questions