How can I code this small piece of code in MEX ?

My code is as follows:

 clc;clear all

%----generate data----%
[X,Y] = meshgrid(-50:2:50, -50:2:50); 
Z = X .* exp(-X.^2 - Y.^2); 
X= X(:);
Y = Y(:);
Z =Z(:);
x = [X Y Z];
%---------------------%

%----- run computation -----%
[m, M] = size (x);
dim = M - 1;

K= zeros (m,m);

for it_dim=1:dim
  tmp = x(:,it_dim+1) * ones(1,m) - ones(m,1) * x(:,it_dim+1)';
  tmp = tmp .* tmp;
  K = K + tmp;
end;
%---------------------------%

can anyone help me out in converting the 'FOR' loop portion of code into MEX code? I am having issues in particular with the following line in conversion to MEX C++:

tmp = x(:,it_dim+1) * ones(1,m) - ones(m,1) * x(:,it_dim+1)';

[James Tursa]

Please use the "{ } Code" button to format your code. Thanks.

Answer 1

[Jan]

Do you want to improve the speed only? Then try this at first:

tmp = bsxfun(@minus, x(:,it_dim+1), x(:,it_dim+1)');

From a mathematical point of view, (a - b)^2 can be written as a^2 - 2*a*b + b^2. This can reduce the number of multiplications massively, because the squaring can be performed for the vectors before inflating them by the multiplication by ONES.

[Kate]

jan, I need to code in Mex since this block of code is part of a larger framework of code which has to be optimized through MEX.

[Jan]

When you speak of optimizing, the main problem here is most likely not the multiplication and summation, but the expensive EXP() function. It does not get clear, if the generation of the test data is a complete dummy or essential part of the problem.

I C you would not replace the vector operations by summing up the single elements.

[Kate]

Jan, the test data component is just a dummy for testing purposes. The 'run computation' component (i.e. the FOR loop) is what I want in MEX form. I am new to MEX programming and having difficulty especially with the indexing,etc in C++. Can you help out?

Answer 2

[Jan]

Perhaps this helps for the conversion to C:

K = zeros(m, m);
for it_dim = 1:dim
   for i1 = 1:m
      K(i1, i1) = K(i1, i1) + (x(i1, it_dim + 1) - x(i1, it_dim + 1)) ^ 2;
      d         = x(i1, it_dim + 1);
      for i2 = i1 + 1:m
         c         = (d - x(i2, it_dim + 1)) ^ 2;
         K(i1, i2) = K(i1, i2) + c;
         K(i2, i1) = K(i2, i1) + c;
      end
   end
end

Btw., the original version is nicer, but this is faster:

 Original/vector operations:  Elapsed time is 0.737606 seconds.
 Elemantary loops, mirroring: Elapsed time is 0.611425 seconds.