Speed up an Enumeration process

Question

After a few days of optimization this is my code for an enumeration process that consist in finding the best combination for every row of W. The algorithm separates the matrix W in one where the elements of W are grather of LimiteInferiore (called W_legali) and one that have only element below the limit (called W_nlegali).

Using some parameters like Media (aka Mean), rho_b_legali The algorithm minimizes the total cost function. In the last part, I find where is the combination with the lowest value of objective function and save it in W_ottimo

As you can see the algorithm is not so "clean" and with very large matrix (142506x3000) is damn slow...So, can somebody help me to speed it up a little bit?

   for i=1:3000
   W = PesoIncertezza * MatriceCombinazioni';
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   W_legali = W;
   W_legali(W<LimiteInferiore) = nan;

   if i==1
        Media = W_legali;
        rho_b_legale = ones(size (W_legali,1),size(MatriceCombinazioni,1));
   else
        Media = (repmat(sum(W_tot_migl,2),1,size(MatriceCombinazioni,1))+W_legali)/(size(W_tot_migl,2)+1);
        rho_b_legale = repmat(((n_b+1)/i),1,size(MatriceCombinazioni,1));
   end

   [W_legali_migl,comb] = min(C_u .* Media .* (1./rho_b_legale) + (1./rho_b_legale) .* c_0 + (c_1./(i * rho_b_legale)),[],2);

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   MatriceCombinazioni_2 = MatriceCombinazioni;
   MatriceCombinazioni_2(sum(MatriceCombinazioni_2,2)<2,:)=[];

   W_nlegali = PesoIncertezza * MatriceCombinazioni_2';
   W_nlegali(W_nlegali>=LimiteInferiore) = nan;

   if i==1
        Media = W_nlegali;
        rho_b_nlegale = zeros(size (W_nlegali,1),size(MatriceCombinazioni_2,1));
   else
        Media = (repmat(sum(W_tot_migl,2),1,size(MatriceCombinazioni_2,1))+W_nlegali)/(size(W_tot_migl,2)+1);
        rho_b_nlegale = repmat(((n_b)/i),1,size(MatriceCombinazioni_2,1));
   end

   [W_nlegali_migliori,comb2] = min(C_u .* Media .* (1./rho_b_nlegale) + (1./rho_b_nlegale) .* c_0 + (c_1./(i * rho_b_nlegale)),[],2);

   z = [W_legali_migl, W_nlegali_migliori];

   [z_ott,comb3] = min(z,[],2);

   %Increasing n_b
   if i==1
       n_b = zeros(size(W,1),1);
   end

   index = find(comb3==1);
   increment = ones(size(index,1),1);
   B = accumarray(index,increment);
   nzIndex = (B ~= 0);
   n_b(nzIndex) = n_b(nzIndex) + B(nzIndex);

   %Using comb3 to find where is the best configuration, is in
   %W_legali or in W_nLegali?

   combinazione = comb.*logical(comb3==1) + comb2.*logical(comb3==2);
   W_ottimo = W(sub2ind(size(W),[1:size(W,1)],combinazione'))';

   W_tot_migl(:,i) = W_ottimo;
   FunzObb(:,i) = z_ott;


   [PesoCestelli] = Simulazione_GenerazioneNumeriCasuali (PianoSperimentale,NumeroCestelli,NumeroEsperimenti,Alfa);
   [PesoIncertezza_2] = Simulazione_GenerazioneIncertezza (NumeroCestelli,NumeroEsperimenti,IncertezzaCella,PesoCestelli);

   PesoIncertezza(MatriceCombinazioni(combinazione,:)~=0) = PesoIncertezza_2(MatriceCombinazioni(combinazione,:)~=0); %updating just the hoppers that has been discharged

end

@chappjc I did it right now and I preallocate `FunzObb` and `W_tot_migl`. I was hoping in something big because right now it's very slow... — gmeroni, Oct 29 '13 at 05:52
@chappjc This is what i get: http://f.cl.ly/items/0i1h1A1A362H1W3N1E22/file20.html (I run it with bad input value for increasing the speed, but normally the `Enumerazione_fo` call is around 5000s or more) — gmeroni, Oct 29 '13 at 07:18

chappjc · Accepted Answer · 2013-10-29T08:13:19.703

1

When you see repmat you should think bsxfun. For example, replace:

Media = (repmat(sum(W_tot_migl,2),1,size(MatriceCombinazioni,1))+W_legali) / ...
    (size(W_tot_migl,2)+1);

with

Media = bsxfun(@plus,sum(W_tot_migl,2),W_legali) / ...
    (size(W_tot_migl,2)+1);

The purpose of bsxfun is to do a virtual "singleton expansion" like repmat, without actually replicating the array into a matrix of the same size as W_legali.

Also note that in the above code, sum(W_tot_migl,2) is computed twice. There are other small optimizations, but changing to bsxfun should give you a good improvement.

The values of 1./rho_b_legale are effectively computed three times. Store this quotient matrix.

edited Oct 29 '13 at 08:13

answered Oct 29 '13 at 07:59

chappjc

30,359
6
75
132

Thanks a lot, `bsxfun` works like a charm! I fixed the `sum(W_tot_migl,2)` and add some code changes. It's a bit faster! If something pop up in your mind feel free to share! Right now the profiler (*the real one*) plot this: http://f.cl.ly/items/3g0M1x0r092h3l0e4643/file18.html – gmeroni Oct 29 '13 at 12:50
One more thing, I tried to use bsxfun. Is this correct? `n_b_legale = repmat((n_b+1),1,size(MatriceCombinazioni,1));` is equal to `n_b_legale = bsxfun(@times,n_b+1,ones(1,size(MatriceCombinazioni,1)));` – gmeroni Oct 30 '13 at 08:24
1

I guess this is a case where you should not think `bsxfun`. :) You don't need to use `bsxfun` with a scalar like `n_b+1`. Instead, do `(n_b+1)*ones(1,size(MatriceCombinazioni,1))`. – chappjc Oct 30 '13 at 15:25

Speed up an Enumeration process

1 Answers1