I have a "super-nested" loop in Matlab to create an array gamma. The loop takes a lot of time to execute. I would like your help to vectorise the loop.

This is the code with explanations:

rng default
clear 

%%%%%%%%%%%%%%%%%Parameters
L=2;
K=3;
n_draws=10^6;
mu_V=zeros(1, L 1);  
Sigma_V=eye(L 1);  
v_draws=mvnrnd(mu_V,Sigma_V,n_draws); %n_drawsx(L 1)
payoff=randn(n_draws, L 1); %n_drawsx(L 1)
v_upp=5;  
v_low=-5; 

%%%%%%%%%%%%%%%%Fill the array gamma
gamma=zeros(K 1,K 1,K 1,L 1,L 1); %allocate space for the array gamma
for k1=1:K 1 
    for k2=1:K 1
        for k3=1:K 1 
            for y=1:L 1  
                for y1=1:L 1
                    tic
                    integr=((nchoosek(K, (k1-1))*(((v_draws(:,1)-v_low).^(k1-1).*(v_upp-v_draws(:,1)).^(K-(k1-1)))./((v_upp-v_low).^K))).*...
                            (nchoosek(K, (k2-1))*(((v_draws(:,2)-v_low).^(k2-1).*(v_upp-v_draws(:,2)).^(K-(k2-1)))./((v_upp-v_low).^K))).*...
                            (nchoosek(K, (k3-1))*(((v_draws(:,3)-v_low).^(k3-1).*(v_upp-v_draws(:,3)).^(K-(k3-1)))./((v_upp-v_low).^K)))).*...
                             mvnpdf(v_draws,mu_V,Sigma_V).*...
                            (payoff(:,y)-payoff(:,y1)); %n_drawsx1  
                    gamma(k1,k2,k3,y,y1)=sum(integr)/n_draws; %average of elements of integr
                    toc
                end
            end
        end
    end
end

For example, with K=3, it takes around 20 sec to execute. In my actual exercise, K=50 and it takes forever.

CodePudding user response：

You can save ~95% runtime from my testing by applying two main types of improvement

1. Compute things in the outermost loop possible, so that equivalent terms are not computed multiple times.
2. Vectorise the inner y1 loop to remove it

All of the deltas between v_upp, v_low, and v_draws are constant because these vectors do not change in the loops, so can all be computed outside of the loops. So can the mvnpdf.

The code then looks something like this:

gamma=zeros(K 1,K 1,K 1,L 1,L 1); %allocate space for the array gamma
v_delta = v_upp-v_low;
v_delta_upp1 = v_upp-v_draws(:,1);
v_delta_upp2 = v_upp-v_draws(:,2);
v_delta_upp3 = v_upp-v_draws(:,3);
v_delta_low1 = v_draws(:,1)-v_low;
v_delta_low2 = v_draws(:,2)-v_low;
v_delta_low3 = v_draws(:,3)-v_low;
MVN = mvnpdf(v_draws,mu_V,Sigma_V);
y1=1:L 1;
for k1=1:K 1 
    term1 = (nchoosek(K, (k1-1))*((v_delta_low1.^(k1-1).*v_delta_upp1.^(K-(k1-1)))./(v_delta.^K)));
    for k2=1:K 1
        term2 = (nchoosek(K, (k2-1))*((v_delta_low2.^(k2-1).*v_delta_upp2.^(K-(k2-1)))./(v_delta.^K)));
        for k3=1:K 1 
            term3 = (nchoosek(K, (k3-1))*((v_delta_low3.^(k3-1).*v_delta_upp3.^(K-(k3-1)))./(v_delta.^K)));
            integr = term1 .* term2 .* term3 .* MVN;
            for y=1:L 1  
                integrp = integr .* (payoff(:,y)-payoff(:,y1)); %n_drawsx1  
                gamma(k1,k2,k3,y,y1)=sum(integrp)/n_draws; %average of elements of integr                
            end
        end
    end
end

You could further improve this by noticing things like

• You're calculating nchoosek of 1:K 1 inside every loop, it might be quicker to calculate an array of the nchoosek values up front and just index into it when looping.
• You're dividing by n_draws and multiplying by the mvnpdf term every loop iteration, is it quicker to just do these after the loop?

You'll get diminishing returns looking at things like this though.