vectorial programming instead of for-loops

Hi All,

I have implemented one function that I would like to improve. I used two for loops and instead of at least one of them I would like to use vectorial programming. (The reason is, that my input xMatrix will contain about 100000 columns and the run time will be very large then.) Unfortunately I don´t have an idea how to use the vectorial programming at all. Could you please give me some advices how can I improve the run time of the function below?

That´s my code:

function f = f_n(xMatrix, MatrixX)
n=size(MatrixX,2);
hn=0.5
for i=1:size(xMatrix,2)
    xi=xMatrix(:,i);
    sume = ones(1, size(xMatrix,2));
    sume(i)=0;
    for j=1:n;
        Xj=MatrixX(:,j);
      g=(xi-Xj)./hn;
      sume(i) = sume(i) + g;
    end
end
f=sume/(n*(hn^3));

Thank you in advance Best Regards Xanka

[Image Analyst]

Why are you preallocating sume with ones(), then assigning some values to it, then completely blowing away all your assignments by setting it equal to an array of all 1's again on the next i'th iteration? You've destroyed all your work! Because of that, I can't really figure out what you want to do. I think conv2() might be able to help but with bad code and no comments, I can't really tell.

Also, could the confusingly named xMatrix and MatrixX be any arbitrary size? Like xMatrix is 5 by 100,000 and MatrixX is 37 by 42? Or do some of their dimensions equal each other? If so, which dimension(s)?

[Cedric]

You should explain what you want to achieve, and give an example, because we cannot figure out based on your code.

[xanka bocem]

Hi all, I will try to explain my problem better. Please find attached also the formula I´m triing to implement.

xMatrix contains 3 rows, 100000 columns and the MatrixX containss 3rows, 200 columns.

What I tried to do is to pick a fixed column (that´s a vector) from the xMatrix and substract from this vector each vector of the MatrixX. The difference is also a vector. Than I should calculate the norm of this calculated vector and sum these norms to a value. This value is the result for the choosen column of the xMatrix. And then the same for the remaining columns of xMatrix. So at the end I should get an array with exact so many values as columns in the xMatrix. That´s why I created at the beginning an array where I overwrite the entries with calculated values. If it isn´t necessary to define such an array I can also omit it.

Here adjusted code (I saw that norm was missing.)

function f = f_n(xMatrix, MatrixX)
n=size(MatrixX,2);
hn=0.5
for i=1:size(xMatrix,2)
    xi=xMatrix(:,i);
    sume = ones(1, size(xMatrix,2));
    sume(i)=0;
    for j=1:n;
        Xj=MatrixX(:,j);
      g= norm((xi-Xj)./hn);
      sume(i) = sume(i) + g;
    end
end
f=sume/(n*(hn^3));

I hope I could explain it better now. Best Regards

Answer 1

[Andrei Bobrov]

hn = .5;
n=size(MatrixX,2);
f = ...
 sum(sqrt(sum((bsxfun(@minus,xMatrix,permute(MatrixX,[1,3,2]))/hn).^2)),3)/(n*(hn^3));

[xanka bocem]

Thank you for the solution! :)

Answer 2

[Cedric]

The first thing to consider is that given the size of matrices that you are dealing with (and the fact that elements are of type double, stored on 8 bytes), storing all components of the differences (in the norm) should take around:

 >> 1e5 * 2e2 * 3 * 8 / 1e6
 ans =
   480

(1e5*2e2 differences between vectors, 3 components per vector, 8 bytes per component), so 480MB of RAM, which is not a lot. It may therefore be more efficient to compute everything simultaneously using arrays with all differences/norms/etc, than to loop over a large dimension (even if you were using smaller arrays this way).

Here is some code where I am using a small test setup that will allow you to follow the computations and check by hand (sorry if there are mistakes, I didn't have much contiguous time). In short, what we do is that we build a giant array of components of differences x-Xi that we use to perform further operations. BSXFUN helps us computing these differences, and we loop over the smallest dimension, which is the number of components: 3. Then we divide by hn, compute norms, sum over the dimension that corresponds to X, and divide by the scalar factor.

 clear
 clc

 % - Small test setup.
 xs = [5, 7; 5, 7; 5, 7] ;
 X  = magic( 3 ) ;
 hn = 0.5 ;

 % - Prealloc 3D array of components; dim 1 associated with vectors of X,
 %   dim 2 with vectors of xs, and dim 3 are for 3 dims/comps per vector.
 n_xs  = size( xs, 2 ) ;
 n_X   = size( X, 2 ) ;
 comps = zeros( n_X, n_xs, 3 ) ;

 % - Build array of differences, looping over dims (3).
 for dimId = 1 : 3
   comps(:,:,dimId) = bsxfun( @minus, xs(dimId,:), X(dimId,:)' )  ;
 end

 % - Divide differences by hn, and compute norms of components 
 %  differences/hn.
 norms = sqrt( sum( (comps/hn).^2, 3 )) ;

 % - Sum over the dimension that corresponds to the number of columns
 %  of X and divide by scalar factor.
 f_xs = 1 / (n_X * hn^3) * sum( norms, 1 ) ;

Running this, we get:

 xs =
     5     7
     5     7
     5     7

 X =
     8     1     6
     3     5     7
     4     9     2

 comps(:,:,1) =
    -3    -1
     4     6
    -1     1
 comps(:,:,2) =
     2     4
     0     2
    -2     0
 comps(:,:,3) =
     1     3
    -4    -2
     3     5

 norms =
    7.4833   10.1980
   11.3137   13.2665
    7.4833   10.1980

 f_xs =
   70.0809   89.7669

So you can check manually that it is correct, or use it to find mistakes.

Hope it helps! Cedric

[Cedric]

PS: I just see Andrei's solution. If you want to understand what happens, I advise you (Xanka) to study my solution with the three 2D singleton expansion (with BSXFUN). But then, Andrei implemented a direct 3D one which makes his solution about 25% faster than mine, so you should implement his solution!

If you understood my simple expansion, Andrei transforms MatrixX with the following:

 permute(MatrixX,[1,3,2])

into an array whose singleton dimension is dim 2 (you would see it along axes x and z in a 3D Cartesian system of axes). Then he makes an expansion applying the subtraction operator with matrix xMatrix, whose singleton dimension is implicitly dim 3. This creates an array of differences which is rotated with respect to mine (the 3 components of what you call vectors are along dim 1 in his case, and along dim 3 in mine). Then he divides by hn and squares, and sums implicitly along dim 1, whereas I do it along dim 3. Etc..

Finally, all this goes quite beyond what people generally call "vectorizing the code" ( ref #1, ref #2 ); so no worry if it looks very complicated to you, because it is not an easy form of vectorization.

[xanka bocem]

Thank you for the detailed explanation! :)