How to vectorize for loops?

Hi Everybody, I have three for loops and their processing is very slow, I need to speed up the process. For that purpose we need to convert it to vectors . Any help will be strongly encouraged. Below is the code:

 for k = 1:size_glcm_3 
    glcm_sum(k) = sum(sum(glcm(:,:,k)));
    glcm(:,:,k) = glcm(:,:,k)./glcm_sum(k); % Normalize each glcm
    glcm_mean(k) = mean2(glcm(:,:,k)); % compute mean after norm
    glcm_var(k)  = (std2(glcm(:,:,k)))^2;

      for i = 1:size_glcm_1  
        for j = 1:size_glcm_2
            p_x(i,k) = p_x(i,k) + glcm(i,j,k); 
            p_y(i,k) = p_y(i,k) + glcm(j,i,k); % taking i for j and j for i
            p_xplusy((i+j)-1,k) = p_xplusy((i+j)-1,k) + glcm(i,j,k);
            p_xminusy((abs(i-j))+1,k) = p_xminusy((abs(i-j))+1,k) +...
            glcm(i,j,k);
        end 
    end
end

All arrays are pre-allocated, size of size_glcm_1 and size_glcm_2 is 512 and size of size_glcm_3 is 1 .

[Jan]

Please edit the question and insert more information. Do not post code lines, which are commented, because this confused the readers. Provide some example data, because the typical dimensions matter, e.g. if size_glcm_1 is 10 and size_glcm_2 is 1e7 or the other way around.

Please run the üprofile to find the bottleneck of the code. If 80% of the processing time are spent inside mean2 optimizing the loops will not help that much.

Are the arrays pre-allocated? Is gclm required as result or is it a temporary variable only? p_y is the tranposed p_x, so do you really need to calculate it?

[azizullah khan]

@Jan Thank you for your response. Yes, i pre-allocated all arrays, both the size of size_glcm_1 and size_glcm_2 is 512 and size_glcm_3 is 1. the time it is taking is due the two internal for loops ( i & j). Waiting for your response. thx

[Image Analyst]

Why not use var() instead of std2() and squaring?

Answer 1

[dpb]

To simplify, I presume copy the plane slice to a 2D array to eliminate the k index from the expressions. The following duplicated your results for a trial array of random values...

px=sum(glcm,2);
py=sum(glcm,1);
N=size_glcm_1-1;
j=0;
for i=-N:N
  j=j+1;
  pxp(j)=sum(diag(fliplr(glcm),i));
end
pxm(1)=sum(diag(glcm));   % is only one main diagonal
for i=1:N
  pxm(i+1)=sum(diag(glcm),-i)+sum(diag(glcm),i); % +/- off-diagonals
end

NB: You may want a temporary for fliplr(glcm); I'm not sure if the JIT optimizer will avoid doing the operation every pass or not; you can test and adjust as seems necessary.

You can also 'spearmint w/ accumarray and friends to see about eliminating the remaining loops; it wasn't patently apparent to me it would help altho for a fixed size you perhaps could build the needed indexing arrays a priori.

ADDENDUM

OK, the accumarray solution isn't as bad as I thought it might have been...see comments below on "how it works".

Alternate solution--

% the preliminaries...
N=size_glcm_1-1;
[i j]=ind2sub(size(glcm),1:numel(glcm));
idx=i+j-1; idx(end)=1;
% calculations
px=sum(glcm,2);
py=sum(glcm,1);
pxp=accumarray(idx.',glcm);
pxm(1)=sum(diag(glcm));   % is only one main diagonal
for i=1:N
  pxm(i+1)=sum(diag(glcm),-i)+sum(diag(glcm),i); % +/- off-diagonals
end

[i j]=ind2sub(size(glcm),1:numel(glcm));
idx=i+j-1; idx(end)=1;
pxp=accumarray(idx.',glcm); % the result

[azizullah khan]

Thank you for your response. But i'm not clear about the above code, Will it perfrom what my code is performing and seconde One, will it decrease the time ??Alsco check the second for loop for j. Please state , i'm waiting, thx

[dpb]

...The following duplicated your results for a trial array..."

What more can I say?

As for the rest, it's your problem, not mine; you'll have to do some timing tests although I'd certainly expect it to outperform the straightahead indexing case altho again TMW has markedly improved performance of loops w/ the JIT optimizer so it's not nearly as certain a foregone conclusion that vectorized code will do better than the straightforward looping solution.

You are aware that your pxm(inus) indexing results in only half the length of the resulting array as does xpp(lus), since the abs() operation folds the two diagonal locations on top of each other? I presume this is intentional?

Oh, one last thing on speed; I'll be quite certain the first two sum expressions will win; whether sum(diag(...)) wins or not is an "exercise for the student" to determine.

Second "oh"--the j was a typo; copied the above loop and pasted/edited in place rather than from the testing at command line and missed fixign it, sorry...

You might find the following at the command line useful for visualizing if you're having trouble following -- use a small array size (like 4x4, say) to illustrate.

>> for i=1:4,for j=1:4,kp=i+j-1;km=abs(i-j)+1;disp([i j kp km]),end,end

That'll be a table of the i,j indices and the "plus/minus" indices for each pair for the array. Consolidating those last two columns shows the subsets that are selected for each term if you've not already done the exercise....

[dpb]

Here's the testing I did at command line to verify what your coding was doing...

% some sample data...
g =[68    66    28    70
    76    18     5    32
    75    71    10    96
    40     4    83     4];
>> px=zeros(1,4);py=px;  % preallocate, compute simple sums
>> for i=1:4,for j=1:4,px(i)=px(i)+g(i,j);py(i)=py(i)+g(j,i);end,end
>> px, py
px =
 232   131   252   131
py =
 259   159   126   202
>> all(px==sum(g,2).')  % check get same results...
ans =
   1
>> all(py==sum(g))
ans =
   1
>>

Both show 'True' so that was easy-peasy...

After the above sample index table was created it was then pretty simple to build the other solution...first your version (remember I just did this at the command line; this is just cut 'n paste from there so the formatting isn't "prettified", only some comments added as went thru the steps what was doing/thinking to solve the problem)--

>> for i=1:4,for j=1:4,
pxp(i+j-1)=pxp(i+j-1)+g(i,j);k=abs(i-j)+1;pxm(k)=pxm(k)+g(j,i);end,end
>> [pxp.'; pxm.']  % answers; note the zeros for the minus case
ans =
  68   142   121   186    46   179     4
 100   397   139   110     0     0     0
>>

After the table noted the diagonals as you've written are the reverse of the normal diagonals; ie, they're +45 not -45 angle. To use diag easily, fliplr g --

>> m=fliplr(g);
>> for i=-3:3,sum(diag(m,i)),end % I elided the other 'ans' text for brevity
ans =
   4
 179
  46
 186
 121
 142
  68
>>

By inspection you observe the same values although in reverse order; if order is important switch the order of the loop from -N:N to N:-1:-N in posted code.

I'll leave the other to you to finish the comparison; the technique should be apparent by now.

[azizullah khan]

Thank you for your response.The rest of data is ok but Now I'm getting error in the last loop.PLease review this section of code.

for i=1:N
pxm(i)=sum(diag(glcm),-i)+sum(diag(glcm),i);
end

I would really appreciate if you can provide proof like the above.You almost solved my problem. thx

[Walter Roberson]

What error is showing up? Please post everything in red.

[dpb]

Oh, ok, I see your problem...in

pxm(i)=sum(diag(glcm),-i)+sum(diag(glcm),i);

You've got the offset argument in the diag call in the wrong place--

pxm(i)=sum(diag(glcm,-i))+sum(diag(glcm,i));

[dpb]

Updated Answer to reflect above oversight and a couple minor typos observed...

[azizullah khan]

    pxm(1)=sum(diag(glcm));   % is only one main diagonal
for i=1:N
  pxm(i+1)=sum(diag(glcm),-i)+sum(diag(glcm),i); % +/- off-diagonals
end
*Error using sum
Dimension argument must be a positive integer scalar within indexing range.*

[dpb]

Same error...see above.

[azizullah khan]

Finally got Luck with replacing the last for loop , but fliplr is time consumable , if any alternative for that would be better in term of computational time.. Final Last Loop:

   for i=1:N,pxm(i+1)=sum(diag(glcm,-i))+sum(diag(glcm,i)),end

[dpb]

"... fliplr is time consumable"

Did you follow the suggestion of making the result of the operation a temporary outside the loop?

Else, I'd have to think about it some; it's the way you've defined the subscripts unless you could live with the alternate diagonals in the "normal" direction.

ADDENDUM

Actually, on reflection, there is a rather cute relationship that can be exploited that should be reasonably fast.

[i j]=ind2sub(size(g),1:numel(g))-1;  % the indices of the array less one
idx=i+j; idx(end)=1;     % the combination of the two, correct last position
pxp=accumarray(idx.',g); % the result

It turns out the "positive" diagonals have corresponding indices whose sum is a constant excepting for the first and last which are combined as the first element; hence the fixup for the last position in the summation.

The index vector can be computed outside any loop, of course, and if they were to always be the same could be precomputed and simply read or passed in saving even the generation time.

One would have to test whether this beats the fliplr or not; I'm guessing as it also eliminates the double call to diag it will...

[azizullah khan]

Thank you @dpb , I didn't test flip outside the loop, No, i tested it outside and result is tremendous, Computataional time is decreased now !!!!!