Matlab taking too much time in computing mex function for large array

Question

I wrote a MATLAB script in which I am passing few scalars and one row vector as input arguments of a mex function and after doing some calculation, it is returning a scalar as output. This process has to be done for all the elements of an array whose size is 1 X 1638400. Below is the corresponding code:

ans=0;
for i=0:1638400-1
    temp = sub_imed(r,i,diff);
    ans  = ans + temp*diff(i+1); 
end

where r,i are scalars, diff is a vector of size 1 X 1638400 and sub_imed is a MEX function which does the below job:

void sub_imed(double r,mwSize base, double* diff, mwSize dim, double* ans)              
{                                                                                           
     mwSize i,k,l,k1,l1;
     double d,g,temp;

     for(i=0; i<dim; i++)
     {   
          k = (base/200) + 1;
          l = (base%200) + 1;
          k1 = (i/200) + 1;
          l1 = (i%200) + 1;

          d = sqrt(pow((k-k1),2) + pow((l-l1),2));

          g=(1/(2*pi*pow(r,2)))*exp(-(pow(d,2))/(2*(pow(r,2))));   

          temp = temp + diff[i]*g;
     }

     *ans  = temp;
}

void mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[]) 
{
    double *diff;           /* Input data vectors */
    double r;               /* Value of r (input) */
    double* ans;            /* Output ImED distance */
    size_t base,ncols;      /* For storing the size of input vector and base */

    /* Checking for proper number of arguments */
    if(nrhs!=3) 
       mexErrMsgTxt("Error..Three inputs required.");

    if(nlhs!=1) 
       mexErrMsgTxt("Error..Only one output required.");

    /* make sure the first input argument(value of r) is scalar */
    if( !mxIsDouble(prhs[0]) || mxIsComplex(prhs[0]) || mxGetNumberOfElements(prhs[0])!=1) 
       mexErrMsgTxt("Error..Value of r must be a scalar."); 

    /* make sure that the input value of base is a scalar */
    if( !mxIsDouble(prhs[1]) || mxIsComplex(prhs[1]) || mxGetNumberOfElements(prhs[1])!=1) 
       mexErrMsgTxt("Error..Value of base must be a scalar."); 

    /* make sure that the input vector diff is of type double */
    if(!mxIsDouble(prhs[2]) || mxIsComplex(prhs[2]))    
       mexErrMsgTxt("Error..Input vector must be of type double.");

    /* check that number of rows in input arguments is 1 */
    if(mxGetM(prhs[2])!=1) 
       mexErrMsgTxt("Error..Inputs must be row vectors."); 

    /* Get the value of r */
    r = mxGetScalar(prhs[0]);
    base = mxGetScalar(prhs[1]);

    /* Getting the input vectors */
    diff = mxGetPr(prhs[2]);
    ncols = mxGetN(prhs[2]);

    /* Creating link for the scalar output */
    plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL);
    ans = mxGetPr(plhs[0]); 

    sub_imed(r,base,diff,(mwSize)ncols,ans);
}

For more details about the problem and the underlining algorithm please follow the thread Euclidean distance between images.

I did a profiling of my MATLAB script and got to know that it is taking 63 sec. just for 387 calls to sub_imed() mex function. So for 1638400 calls to sub_imed, ideally it will take around 74 hours which is just too long.

Can someone please help me to optimize the code by suggesting some alternative ways to reduce the time taken for computation.

Thanks in advance.

Answer 1

I ported your code back to MATLAB and made some small adjustements, while the results should stay the same. I introduced the following constants:

N = 8192;
step = 0.005;

Note that N / step = 1638400. With that, you can rewrite your variable k (and rename it to baseDiv):

baseDiv = 1 + (0 : step : (N-step)).';

i.e. it is 1:8193 in steps of 0.005. Similarly, l is 1:200 (=1:(1/0.005)), repeated 8192 times in a row, which is (now called baseMod):

baseMod = (repmat(1:1:(1/step), 1, N)).';

Your variables k1 and l1 are simply the ith element of k and l, i.e. baseDiv(i) and baseMod(i).

With the vectors baseDiv, and baseMod, one can calculate d, g and your temporary variable tmp with

d = sqrt((baseDiv(k)-baseDiv).^2 + (baseMod(k)-baseMod).^2);
g = 1/(2*pi*r^2) * exp(-(d.^2) / (2*r^2));
tmp = sum(diffVec .* g);

We can put this into your MATLAB for loop, so the whole program becomes

% Constants
N = 8192;
step = 0.005;

% Some example data
r = 2;
diffVec = rand(N/step,1);

base = (0:(numel(diffVec)-1)).';    
baseDiv = (1:step:1+N-step).';
baseMod = (repmat(1:1:(1/step), 1, N)).';

res = 0;
for k=1:(N/step)
    d = sqrt((baseDiv(k)-baseDiv).^2 + (baseMod(k)-baseMod).^2);
    g = 1/(2*pi*r^2) * exp(-(d.^2) / (2*r^2));
    tmp = sum(diffVec .* g);
    res = res + tmp * diffVec(k);
end

By eliminating the inner for loop and calculating it in a vectorized fashion, you still need 11 sec for 1000 iterations, resulting in a total runtime of 5 hours. Still - a speed-up of more than 10x. To get an even higher speed-up, you have two possibilities:

1) Complete vectorization: You can easily vectorize the remaining for-loop by using bsxfun(@minus, baseDiv, baseDiv.') and suming over the columns to calculate all values at the same time. Unfortunately we have a small problem: a 1638400-by-1638400 double matrix would take up 20TB of RAM, which - I assume - you don't have in your Laptop ;-)

2) Less samples: You are doing some mathematical transform with a resolution of step=0.005. Check if you really, really need this precision! If you take 1/10 of the precision: step=0.05, you are 100 times faster, and are finished within 3 minutes!

Answer 2

• Did you compile with the optimisation flag? (See flag -O)
• Instead of using pow(x,2) in C/C++ you should write x*x
• I am almost certain that the reason why it is slow is not because of the code you've shown, but because of what you do in the mexFunction. If I had to guess I'd say you're pointlessly copying the memory in diff, but we need to see the entire Mex function to be sure.

---

Try the following C++ code using mex -O myfile.cpp:

void sub_imed( double r, size_t base, const double *diff, size_t dim, double& ans)
{
    double d, g;

    // these need to be double to avoid underflow
    double k = base / 200;
    double l = base % 200;

    r = 2*r*r;
    for(; dim; --dim, ++diff )
    {
        d = k - i/200;
        g = l - i%200;

        ans += (*diff) * exp( - (d*d + g*g)/r ) / (pi*r);
    }
}

void mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[])
{
    /* Checking for proper number of arguments */
    if(nrhs!=3)
        mexErrMsgTxt("Error..Three inputs required.");

    if(nlhs!=1)
        mexErrMsgTxt("Error..Only one output required.");

    /* make sure the first input argument(value of r) is scalar */
    if( !mxIsDouble(prhs[0]) || mxIsComplex(prhs[0]) || mxGetNumberOfElements(prhs[0])!=1 )
        mexErrMsgTxt("Error..Value of r must be a scalar.");

    /* make sure that the input value of base is a scalar */
    if( !mxIsDouble(prhs[1]) || mxIsComplex(prhs[1]) || mxGetNumberOfElements(prhs[1])!=1 )
        mexErrMsgTxt("Error..Value of base must be a scalar.");

    /* make sure that the input vector diff is of type double */
    if( !mxIsDouble(prhs[2]) || mxIsComplex(prhs[2]) )
        mexErrMsgTxt("Error..Input vector must be of type double.");

    /* check that number of rows in input arguments is 1 */
    if( mxGetM(prhs[2])!=1 )
        mexErrMsgTxt("Error..Inputs must be row vectors.");

    /* Get the value of r */
    double r    = mxGetScalar(prhs[0]);
    size_t base = static_cast<size_t>(mxGetScalar(prhs[1]);

    /* Getting the input vectors */
    const double *diff = mxGetPr(prhs[2]);
    size_t nrows = static_cast<size_t>(mxGetN(prhs[2]));

    /* Creating link for the scalar output */
    plhs[0] = mxCreateDoubleScalar(0.0);
    sub_imed( r, base, diff, nrows, *mxGetPr(plhs[0]) );
}