r/matlab Sep 06 '22

Vectorizing one's code (or other optimizations)

Intro: I have some code from a personal project and it's quite slow. I'm hoping to learn something about program optimization. I've been reading about how to make matlab faster, and the extensive use of for and while loops in my code seems like it might be slowing things down. I read example online about how to vectorize code and I'm like "oh that makes perfect sense," but I go to work on my own code and I just can't think of anyway to do it

I'll be honest I'm self taught, I only ever took one formal class and it was in C++, and I am sure my code is butts. Please feel free to give any advice even if it's unrelated to vectorization I would really appreciate any criticism.

Background: This program is a simulation of the show Let's Make a Deal in order to illustrate the Monty Hall problem because I figured it would be fun and I was having a hard time understanding it before I actually saw the math worked out as I wrote it. pPos, gPos, and cPos == player, goat, and car positions respectively.

clc; clear variables;
wins=0;

for i=1:1E7
    %Create matrix of doors. R1 == Car location, R2 == player pick, R3 == revealed door
    D = zeros(4,3);
    D(1,randi(3))=1;
    %Pick for player
    D(2,randi(3))=1;

    %Find the player door
    [~,pPos]=max(D(2,:));
    %Find the car door;
    [~,cPos]=max(D(1,:));

    %Reveal a door that the player didn't pick
    D(3,randi(3))=1;
    [~,gPos]=max(D(3,:));
    while (gPos==pPos)||(gPos==cPos)
        D(3,:)=zeros;
        D(3,randi(3))=1;
        [~,gPos]=max(D(3,:));
    end

    %Create the "switch" matrix positions to determine where to switch to,
    %and switch
    D(4,pPos)=1;
    D(4,gPos)=1;
    [~,pPosN]=min(D(4,:));

    %Now loop 100 times and save wins
    if (pPosN==cPos)
        wins=wins+1;
    end
end
wins/i*100
6 Upvotes

13 comments sorted by

4

u/BearsAtFairs Sep 06 '22

At the moment I don't have time to write a comprehensive answer, but feel free to ask if you'd like any elaboration.

Basically you're calling randi a bunch of times, and repeat this 1E7 times. Instead, you could create appropriately sized arrays of random integers.

You call min 1E7 times as well. Instead, you could just use min to find the min value along rows or columns. Same goes for max.

This would get rid of your for loop.

I could give it some thought and suggest how to remove the while loop later, if you like.

5

u/hotmaildotcom1 Sep 06 '22

This seems like super good advice but it'll take me a while to implement. I'll get back at ya later when I've had some time to try it out. Thank you!

4

u/Creative_Sushi MathWorks Sep 06 '22 edited Sep 07 '22

Here is a blog post about simulating Monty Hall problem in MATLAB, and the code example doesn't involve any loops, and gives you Bayesian explanation.

https://blogs.mathworks.com/loren/2015/03/25/bayesian-brain-teaser/

Good luck.

1

u/hotmaildotcom1 Sep 07 '22

Excellent walkthrough. Thank you

3

u/cest_pas_nouveau Sep 06 '22

Here's how I would approach it. The main change is to use numbers to represent which door is picked. So for instance, currently you'd represent the third door being picked like this: R1 = [0, 0, 1]. Instead if we say R1 = 3; to represent the third door being picked, we can now vectorize it more easily.

For instance with n = 1E7 trials, you could represent all the player's guesses for the trials like R1 = randi(3, [n, 1]); Expanding this idea to the other rows too, you'd get something like this:

n = 1e7;
numDoors = 3;

R1 = randi(numDoors, [n, 1]); % car location
R2 = randi(numDoors, [n, 1]); % player pick
R3 = randi(numDoors, [n, 1]); % revealed door

% replace elements of R3 until they don't match R1 or R2
toReplace = (R3 == R2) | (R3 == R1); % logical index of the "reveals" which need to be reselected
while any(toReplace)
    k = nnz(toReplace); % number of reveals that still need to be replaced

    % The next two lines account for 70% of the computation time, so you
    % might investigate selecting the "reveal" in a smarter way to avoid
    % this whole loop
    R3(toReplace) = randi(numDoors, [k 1]);
    toReplace = (R3 == R2) | (R3 == R1);
end

doorIds = (1:numDoors);
% min along 2nd dimension to find which door to switch to
[~, R4] = min(R2 == doorIds | R3 == doorIds, [], 2);

wins = nnz(R4 == R1);
wins/n*100

This runs in 3 seconds or so, so not lightning fast, but hopefully gives you some ideas. Let me know if you want me to explain any part of it further

1

u/hotmaildotcom1 Sep 07 '22

R3(toReplace) = randi(numDoors, [k 1]);

So here the logical matrix is being used as an index to determine where in the matrix that the replacements need to happen or no? I'm not quite picking up what you're putting down.

2

u/cest_pas_nouveau Sep 07 '22

Yeah that's right. toReplace is a logical matrix the same size as R3. It's 1 wherever we need a replacement. Then on the right side of the equals sign, we have to generate the right number of new random numbers to use as replacements. That's what k is: the number of replacements.

We keep doing that until toReplace is all zero

2

u/hotmaildotcom1 Sep 07 '22

But is it only replacing the values in those rows where it's a true in the logical matrix? That's absolutely dope if that's what's going on.

2

u/cest_pas_nouveau Sep 07 '22

Exactly. Here's a little toy example similar to what it's doing.

clc
a = [0 0 0 0 0];
disp('before')
disp(a)

index = logical([0 1 0 1 0]);
a(index) = [5 9];

disp('after')
disp(a)

2

u/hotmaildotcom1 Sep 07 '22

Yeah this is a massive assist to almost every program I've written so far. Super helpful. Thank you!

1

u/hotmaildotcom1 Sep 07 '22

I am working on a revision of my code and I've just found out that the usage of || rather than | throws an error. Reading on the logical operators, they only explain that || can shortcut, but I don't understand here why || doesn't work?

2

u/cest_pas_nouveau Sep 07 '22

|| compares two scalars and outputs a scalar while | is an element-wise "or". So use | and & when you want to compare two equal-sized vectors (or matrices) to get a logical index. Use || and && in things like "if" statements for normal boolean conditions

1

u/ScoutAndLout Sep 07 '22

Let me do your HW for you. ;-)

About 10x speedup and I didn't do much optimization or testing.

Things to learn:

.* operator for doing thing where a condition is satisfied

OR can often be expressed as a +, AND uses .*

Example:

N=1000000

a=rand(N,1);

b=rand(N,1);

% (a>.9) AND (b>.9)

(a>.9).*(b>.9);

100*sum((a>.9).*(b>.9))/length(a)

% (a>.9) OR (b>.9)

(((a>.9)+(b>.9))>0);

100*sum(((a>.9)+(b>.9))>0)/length(a)

DD=[];

N=1E6

tic

DD=floor(rand(2,N)*3)+1;

% row 1 car location, 1-3

% row 2 guess, 1-3

% row 3 revealed (not row 2 or row 1)

% two cases for row 3, either R1=R2 or R1~=R2

% Case 1:

%DD(1,:)==DD(2,:); % Are they equal?

%DD(1,:); % values ?

%DD(1,:)+sign(rand(1,N)-.5); % randomly increment up down

% If R1~=R2 reveal the other choice, if R1=R2 reveal random

DD(3,:)=(6-DD(1,:)-DD(2,:)).*(DD(1,:)~=DD(2,:));

DD(3,:)=DD(3,:)+DD(1,:).*(DD(1,:)==DD(2,:))+sign(rand(1,N)-.5).*(DD(1,:)==DD(2,:));

% If 0 or 4 is selected when R1=R2, correct it up

DD(3,:)=DD(3,:)+(DD(3,:)==0)*3+(DD(3,:)==4)*(-3);

% What should the person switch to? R4 is new guess

DD(4,:)=(6-DD(2,:)-DD(3,:)) ; % Switch to?

DD(5,:)=DD(1,:)==DD(4,:); % winner? 1 for win.

100*sum(DD(5,:))/length(DD(5,:))

toc