current community

your communities

more stack exchange communities

Stack Exchange
tour help
Take the 2-minute tour ×
Code Review Stack Exchange is a question and answer site for peer programmer code reviews. It's 100% free, no registration required.
up vote 1 down vote favorite

I have a large matrix samples, I want to make another matrix sampleProb which is the same size, but modified as below: 

%samples is a 1000*1000 matrix, M = 1000
sampleProb = zeros(1000);
for i = 1:size(samples, 2)
  for j = 1:M
    for k = 1:M
      if (k ~= j)
        sampleProb(j, i) = sampleProb(j, i) + normcdf(samples(j, i) - samples(k, i));
      end
    end
  end
end
sampleProb = sampleProb./M;

I know it's best to avoid for loops when possible in Matlab and operations should be vectorised, but I can't figure out what to do. How can I optimise this for performance?

share|improve this question
    
Please retitle your post based on the task that your code performs, rather than how you wish to rework it. More background about what you are trying to accomplish would also be appreciated. –  200_success♦ Feb 22 at 17:31

1 Answer 1

up vote 1 down vote accepted

Introduction and solution code

Since you mentioned making the output matrix sampleProb, I am assuming it to be initialized as all zeros. So, here's the vectorized implementation with the mighty bsxfun -

%// Get the subtractions
sub1 = bsxfun(@minus,samples,permute(samples,[3 2 1]))

%// Calculate the sum of normcdf's in a vectorized fashion
sampleProb = bsxfun(@minus,sum(normcdf(sub1),3),normcdf(sub1(1,:,1)))./M

Benchmarking

Benchmarking Code with M = 100 -

M = 100;
samples =  rand(M,M);
sampleProb1 =  zeros(M,M);

disp('---------------------------------------- With Original Approach')
tic
for i = 1:size(samples, 2)
    for j = 1:M
        for k = 1:M
            if (k ~= j)
                sampleProb1(j, i) = sampleProb1(j, i) + normcdf(samples(j, i) - samples(k, i));
            end
        end
    end
end
sampleProb1 = sampleProb1./M;
toc, clear sampleProb1 i j k

disp('---------------------------------------- With Proposed Approach')
tic
%// Get the subtractions
sub1 = bsxfun(@minus,samples,permute(samples,[3 2 1]));

%// Calculate the sum of normcdf's in a vectorized fashion
sampleProb = bsxfun(@minus,sum(normcdf(sub1),3),normcdf(sub1(1,:,1)))./M;
toc

Runtime results -

---------------------------------------- With Original Approach
Elapsed time is 21.425729 seconds.
---------------------------------------- With Proposed Approach
Elapsed time is 0.046284 seconds.
share|improve this answer
    
Wow, that is fantastic, I did not expect that much of an improvement. Thanks! –  texasflood Mar 3 at 22:19
    
The output from both methods produces exactly the same result, save for a tiny difference which is just insignificant quantisation noise (max(max(sampleProb1 - sampleProb))=2.2204e-16) , very impressive! And yes, you were right, sampleProb is initialised to zero, I'll add that to the question –  texasflood Mar 3 at 22:25

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.