Optimize this matrix process

I wrote this function that has some heavy matrix multiplications. I'm wondering if something can be done in order to optimize it since it's take too long to compute. This function optimized by GA toolbox. (`position_matrix` is usually `[200 x 8]`)

PROFILER OUTPUT: http://cl.ly/0l2t0r1t2Z0L/Profiler.pdf

function [f] = MHW(position_matrix)

global W_best_tot
global Function

nComb = size(position_matrix,1);
dim = size(position_matrix,2);

CombinantionMatrix = combn([0 1],dim);
H_e = 2;
CombinantionMatrix(sum(CombinantionMatrix,2)<H_e,:)=[];
WT = 250;
WL = 241;
nSim = 3000;
CellUncertainty = 0.5 * 10^(-6);

%// Define Product Hopper allocation
P = 3;                                          %// number of products
HA = 4; WT_A = 130;
HB = 2; WT_B = 30;
HC = 2; WT_C = 90;

H = HA + HB + HC;

% OBJECTIVE FUNCTION COST COEFFICIENTS
Cu_A = 0.03; c_p = 0.6; c_f = 734; c_l = 3.2;
Cu_B = 0.09;
Cu_C = 0.04;

W_best_tot = zeros(nComb,nSim);

HopperWeight_A = position_matrix(:,1:HA);
HopperWeight_B = position_matrix(:,HA+1:HA+HB);
HopperWeight_C = position_matrix(:,HA+HB+1:HA+HB+HC);
CombinantionMatrix_A = CombinantionMatrix(:,1:HA);
CombinantionMatrix_B = CombinantionMatrix(:,HA+1:HA+HB);
CombinantionMatrix_C = CombinantionMatrix(:,HA+HB+1:HA+HB+HC);

for ii=1:nSim

      W_A = HopperWeight_A * CombinantionMatrix_A.';
      W_B = HopperWeight_B * CombinantionMatrix_B.';
      W_C = HopperWeight_C * CombinantionMatrix_C.';

      %ADD ERROR
      RandomizeUncertainty = randn(nComb, dim) * CellUncertainty;
      UncertainWeight = abs((1+RandomizeUncertainty) .* position_matrix);
      UncertainWeight = UncertainWeight * CombinantionMatrix.';

      UncertainWeight_A = UncertainWeight(:,1:HA);
      UncertainWeight_B = UncertainWeight(:,HA+1:HA+HB);
      UncertainWeight_C = UncertainWeight(:,HA+HB+1:HA+HB+HC);

      W_Unc_A = UncertainWeight_A * CombinantionMatrix_A.';
      W_Unc_B = UncertainWeight_B * CombinantionMatrix_B.';
      W_Unc_C = UncertainWeight_C * CombinantionMatrix_C.';

      [~,comb_A] = min(abs(W_Unc_A - WT_A),[],2);
      [~,comb_B] = min(abs(W_Unc_B - WT_B),[],2);
      [~,comb_C] = min(abs(W_Unc_C - WT_C),[],2);

      W_best_A = W_A(sub2ind(size(W_A),1:size(W_A,1),comb_A.')).';
      W_best_B = W_B(sub2ind(size(W_B),1:size(W_B,1),comb_B.')).';
      W_best_C = W_C(sub2ind(size(W_C),1:size(W_C,1),comb_C.')).';

      W_best_ABC(:,:,ii) = [W_best_A W_best_B W_best_C];
      W_best_tot(:,ii) = sum(W_best_ABC(:,:,ii),2);

end

n_b = sum(logical(W_best_tot >= WL),2);
n_bA = sum(logical(reshape(W_best_ABC(:,1,:), nComb, nSim) >= WT_A),2);
n_bB = sum(logical(reshape(W_best_ABC(:,2,:), nComb, nSim) >= WT_B),2);
n_bC = sum(logical(reshape(W_best_ABC(:,3,:), nComb, nSim) >= WT_C),2);

n_b = logical(W_best_tot >= WL);
n_bA = logical(reshape(W_best_ABC(:,1,:), nComb, nSim) >= WT_A);
n_bB = logical(reshape(W_best_ABC(:,2,:), nComb, nSim) >= WT_B);
n_bC = logical(reshape(W_best_ABC(:,3,:), nComb, nSim) >= WT_C);

n_b_global = sum([n_b .* n_bA .* n_bB .* n_bC],2);

Avg   = mean(W_best_tot,2);
StdDev  = sqrt(var(W_best_tot.')).';

GAvg    = sum(W_best_tot.*(W_best_tot > WL),2)./sum(W_best_tot > WL,2);
GAvg_A  = sum(reshape(W_best_ABC(:,1,:), nComb, nSim).*(reshape(...
            W_best_ABC(:,1,:), nComb, nSim) > WT_A),2)./sum(reshape(...
            W_best_ABC(:,1,:), nComb, nSim) > WT_A,2);
GAvg_B  = sum(reshape(W_best_ABC(:,2,:), nComb, nSim).*(reshape(...
            W_best_ABC(:,2,:), nComb, nSim) > WT_B),2)./sum(reshape(...
            W_best_ABC(:,2,:), nComb, nSim) > WT_B,2);
GAvg_C  = sum(reshape(W_best_ABC(:,3,:), nComb, nSim).*(reshape(...
            W_best_ABC(:,3,:), nComb, nSim) > WT_C),2)./sum(reshape(...
            W_best_ABC(:,3,:), nComb, nSim) > WT_C,2);

f = Cu_A .* GAvg_A + Cu_B .* GAvg_B + Cu_B .* GAvg_B + (nSim./...
        n_b_global) .* c_p + (c_f./n_b_global) + ((nSim - n_b_global)./...
        n_b_global) .* c_l;

Function = f;

[dpb]

The first rule of opimization is "profile, profile, profile"...until then you really don't know where, specifically to look if there is a particular bottleneck.

I don't know if the JIT is smart enough to recognize it, you might try in such expressions as

HopperWeight_A(CombinantionMatrix_A(comb_A,:)~=0) = ...
        HopperWeight_2_A(CombinantionMatrix_A(comb_A,:)~=0);

to take the logical outside just in case it does actually get done twice...

idx=CombinantionMatrix_A(comb_A,:)~=0);
HopperWeight_A(idx) = HopperWeight_2_A(idx);

I really don't know; just a guess that shouldn't hurt, might help...

[Gianmarco]

I've updated the code! This is the profiler output: http://cl.ly/0l2t0r1t2Z0L/Profiler.pdf

[dpb]

So, did the suggested have any effect?

The following lines are roughly equivalent cost and about 40% of total altho the following includes some others besides the top 5 that are essentially equivalent so one presumes they'll add a similar relative fraction making the total of this section closer to 70% overall--

W_A = HopperWeight_A * CombinantionMatrix_A.';
W_B = HopperWeight_B * CombinantionMatrix_B.';
W_C = HopperWeight_C * CombinantionMatrix_C.';
...
UncertainWeight = UncertainWeight * CombinantionMatrix.';
...
W_Unc_A = UncertainWeight_A * CombinantionMatrix_A.';
W_Unc_B = UncertainWeight_B * CombinantionMatrix_B.';
W_Unc_C = UncertainWeight_C * CombinantionMatrix_C.';

The other thing that strikes me that I also don't know about is that you're doing the transpose every time inside this deep loop. Can you rearrange storage to avoid that (or at least run an outside timing run to see if again the JIT is smart enough to be able to get around the actual work involved)?

[Gianmarco]

Some improvement were made! :) These are the new profiler files: This is the Main Profiler and these are: the MHW function and the Weight function . The main problems are in the Weight since it's called 3000times for each GA generation...Above the Weight function code:

function [ HopperWeight UncertainWeight ] = Weight( nComb, H, alpha, AllComb, CellUncertainty )

      HopperWeight = randn(nComb,H) * alpha;
      HopperWeight = (1+HopperWeight) .* AllComb;
      idx = HopperWeight<0;
      random = randn(sum(idx(:)), 1);
      HopperWeight(idx) = random .* ((1+alpha) .* AllComb(idx));

      %ADD ERROR
      RandomizeUncertainty = randn(nComb,H) * CellUncertainty;
      UncertainWeight = abs((1+RandomizeUncertainty) .* HopperWeight);
      end

The new MHW function can be found here