MATLAB: Optimize/speed up a big and slow matrix operation with addition and bsxfun

bsxfunelement wise operationsexecution-timematrix additionoptimizationperformance

Hi. I have a fairly long line of code in my script that takes about 7 seconds to run, it was originally three separate calculations but are now done in one line that looks like this:

X = bsxfun(@times,reshape(bsxfun(@times,A,B),[1440,1,2251]),(1-C)./2)...
    +bsxfun(@times,reshape(E,1440,1,2251),D)...
    +bsxfun(@times,reshape(F,1440,1,[]),((1+C)/2));

Since I need to run it several tens of thousands times it is causing a fair amount of anxiety in my life at the moment. I don’t know how I would go about speeding it up further or if it is even possible. So I would really appreciate it if any of you optimization gurus here could give me some advice on how and if I can go about making this faster.

The input variables are of size:

A = 2251x1
B = 1440x2251
C = 1x181
D = 1440x181
E = 1440x2251
F = 1440x2251

Thanks!

Best Answer

Assuming you have a C compiler available, here is a simple C mex routine to do the calculation. Could possibly be sped up even more with multi-threading (e.g., OpenMP) as well but I did not include that code below. Let me know if you have an OpenMP compilant compiler since that addition would be fairly simple.

As is (with no C-mex multi-threading) here is an example run on my machine:

Elapsed time is 6.655025 seconds. % M-code with bsxfun
Elapsed time is 1.700077 seconds. % C-mex code (not multi-threaded)

CAUTION: Code below is bare bones with no argument checking.

// X =   bsxfun(@times,reshape(bsxfun(@times,A,B),[1440,1,2251]),(1-C)./2)...
//     + bsxfun(@times,reshape(E,1440,1,2251),D)...
//     + bsxfun(@times,reshape(F,1440,1,[]),((1+C)/2));
// 
// The input variables are of size:
// 
// A =    1 x       2251
// B = 1440 x       2251
// C =    1 x 181
// D = 1440 x 181
// E = 1440 x       2251
// F = 1440 x       2251
//
// Output:
//
// X = 1440 x 181 x 2251
//
// Calling sequence:
//
// X = this_file_name(A,B,C,D,E,F);
#include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    mwSize i, j, k, m, n, p, km, ikm, jm;
    mwSize dims[3];
    double *A, *B, *C, *D, *E, *F, *X, *Cm, *Cp;
    double a, cp, cm;
      m = mxGetM(prhs[3]);
      n = mxGetN(prhs[3]);
      p = mxGetN(prhs[4]);
      A = mxGetPr(prhs[0]);
      B = mxGetPr(prhs[1]);
      C = mxGetPr(prhs[2]);
      D = mxGetPr(prhs[3]);
      E = mxGetPr(prhs[4]);
      F = mxGetPr(prhs[5]);
      dims[0] = m;
      dims[1] = n;
      dims[2] = p;
      plhs[0] = mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL);
      X = mxGetPr(plhs[0]);
      Cm = (double *) mxMalloc(n*sizeof(*Cm));
      Cp = (double *) mxMalloc(n*sizeof(*Cp));
      for( j=0; j<n; j++ ) {
          Cm[j] = (1.0 - C[j]) / 2.0;
          Cp[j] = (1.0 + C[j]) / 2.0;
      }
      for( k=0; k<p; k++ ) {
          a = A[k];
          km = k * m;
          for( j=0; j<n; j++ ) {
              jm = j * m;
              cm = Cm[j] * a;
              cp = Cp[j];
              for( i=0; i<m; i++ ) {
                  ikm = i+km;
                  *X++ = B[ikm] * cm + E[ikm] * D[i+jm] + F[ikm] * cp;
              }
          }
      }
      mxFree(Cp);
      mxFree(Cm);
}

Related Solutions

MATLAB: Can A = A + B’*B be sped up somehow? It is seriously bottlenecking the for-loop

If you have a C compiler available, another option is to call a BLAS routine from within a mex function to do this calculation. It may not save operation time on the B'*B part (not sure if MATLAB is smart enough to call DSYRK instead of DGEMM), but it may save on the time spent in the other parts. That is because A = A + B'*B can be done with one BLAS routine call (this type of update is built into the BLAS library code), as opposed to doing B'*B first into temporary storage, then adding it to A, then free'ing the temporary storage. E.g., the doc for DSYRK shows that it can do the following operation in one call, without forming B' explicitly first:

A := alpha * B^T * B + beta * A

So for your case, alpha = 1 and beta = 1. Also, by utilizing this BLAS call everything can be done in-place (i.e., no need to temporarily form the right-hand-side result first and then assign it to A).

Let me know if you want to explore this route and if so then I can code it up for you.

EDIT: 9/4/2015

SORRY this has taken so long! I wrote the code actually fairly quickly but had problems in debugging. Code worked in 32-bit but not in 64-bit, but identical code I had in another program worked in both 32-bit and 64-bit. Finally I figured out it was because for this new program I was forgetting to include the -largeArrayDims compilation flag in the mex command (needed to get the mwSignedIndex macro to return the correct type). In any event, here is the code (I changed the integers to ptrdiff_t to match blas.h). Hopefully it will shave some execution time off of your runs.

Two things to keep in mind. First, it does all the calculations IN-PLACE! So make sure your accumulator matrix ("A" in your example) is not shared. If you set it with the zeros function prior to the loop, then it will not be shared and you will be OK. Second, it only does one half of the matrix calculations when using two arguments. After you are done with the loop, then call it once with only one argument to fill in the other half. E.g., A calling sequence to accomplish your above example would be:

A = zeros(6588,6588);
for i =1:1000
    % B changes to a new value here and is a 27x6588 matrix
    dsyrk(A,B); % Does A = A + B'*B in-place, but lower triangle only
end
dsyrk(A); % Fills in the upper triangle from the lower triangle

And here is the mex routine dsyrk.c

/*--------------------------------------------------------------------------------------
 * DSYRK does the operation C = C + A' * A in place
 *
 * Syntax:  dsyrk(C,A) --> does C = C + A' * A (only lower triangle part)
 *          dsyrk(C)   --> fills the upper triangle with the lower triangle
 *
 * C = a real double full N x N matrix
 * A = a real double full K x N matrix
 *
 * The intent is to first initialize C, then call dsyrk(C,A) in a loop for various A
 * matrices, then after the loop call dsyrk(C) once to fill in the upper triangle part.
 * The code does the operations on C in place, so it is up to the user to make sure
 * that C is not shared with any other variable prior to calling dsyrk.
 *
 * The code uses ptrdiff_t for the integer types because that is what the header file
 * blas.h uses for the dsyrk function arguments.
 *
 * Programmer: James Tursa
 *-------------------------------------------------------------------------------------- */
#include "mex.h"
#include "blas.h"
void xFILLPOS(double *Cpr, ptrdiff_t n);
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    double alpha = 1.0, beta = 1.0;
    double *Apr, *Cpr;
    char uplo = 'L';
    char trans = 'T';
    ptrdiff_t m, n, k, p, lda, ldc;
      if( nlhs ) {
          mexErrMsgTxt("Too many outputs ... this routine does all operations in-place");
      }
      if( nrhs == 2 ) {
          Cpr = mxGetPr(prhs[0]);
          m = mxGetM(prhs[0]);
          n = mxGetN(prhs[0]);
          Apr = mxGetPr(prhs[1]);
          k = mxGetM(prhs[1]);
          p = mxGetN(prhs[1]);
          if( m != n || !mxIsDouble(prhs[0]) || mxIsComplex(prhs[0]) || mxIsSparse(prhs[0]) ||
              mxGetNumberOfDimensions(prhs[0]) != 2 ) {
              mexErrMsgTxt("1st input must be real double non-sparse square matrix");
          }
          if( p != n || !mxIsDouble(prhs[1]) || mxIsComplex(prhs[1]) || mxIsSparse(prhs[1]) ||
              mxGetNumberOfDimensions(prhs[1]) != 2 ) {
              mexErrMsgTxt("2nd input must be real double non-sparse matrix compatible with 1st input");
          }
          lda = k;
          ldc = n;
          dsyrk( &uplo, &trans, &n, &k, &alpha, Apr, &lda, &beta, Cpr, &ldc );
      } else if( nrhs == 1 ) {
          Cpr = mxGetPr(prhs[0]);
          m = mxGetM(prhs[0]);
          n = mxGetN(prhs[0]);
          if( m != n || !mxIsDouble(prhs[0]) || mxIsComplex(prhs[0]) || mxIsSparse(prhs[0]) ||
              mxGetNumberOfDimensions(prhs[0]) != 2 ) {
              mexErrMsgTxt("1st input must be real double non-sparse square matrix");
          }
          xFILLPOS(Cpr,n);
      } else {
          mexErrMsgTxt("Need exactly 1 or 2 inputs");
      }
  }
/*--------------------------------------------------------------------------------------
 * Fill the upper triangle with contents of the lower triangle
 *-------------------------------------------------------------------------------------- */
void xFILLPOS(double *Cpr, ptrdiff_t n)
{
    double *source, *target;
    register ptrdiff_t i, j;
      source = Cpr + 1;
      target = Cpr + n;
      for( i=1; i<n; i++ ) {
          for( j=i; j<n; j++ ) {
              *target = *source;
              target += n;
              source++;
          }
          source += i + 1;
          target = source + n - 1;
      }
  }

MATLAB: How to shrink a pre-allocated array in a C-Mex

I am unaware of specific documentation on this, but Part 1 should be sufficient (and safe) as long as you don't care about the wasted memory at the end hanging around (and this is extremely fast of course). I have done this many times myself and have never run into an issue. I have even done this with only one copy of a group of shared data copies (different mxArray variables having different number of elements but pointing to the same data block). I do not believe the MATLAB Memory Manager uses the size information in a mxArray for anything. After all, it is perfectly legal to do this:

mxArray *mx;
double *pr;
mx = mxCreateDoubleMatrix( 0, 0, mxREAL );
pr = mxMalloc( 1000 * sizeof(*pr) );
mxSetPr(mx,pr);
mxSetM(mx,1);
mxSetN(mx,2);
// do stuff
mxSetN(mx,20);
// do more stuff
// etc.

So here you also have a mismatch between the allocated memory and the mxArray size information, but everything is OK as long as the memory is large enough to encompass the sizes present.

What I have had issues with in the past is with mxRealloc not always freeing up the memory when going from a large block to a small block. Or, at least the memory didn't always show up as available after the mxRealloc call (maybe the memory manager has it somehow marked internally as available but I didn't have access to the details so I couldn't really tell).

Side Issue:

The following comparison is technically non-conforming:

q != p

If p was free'd by the mxRealloc call, you can't legally use it in downstream code, even for comparison purposes (although you can probably get away with it on the machines we typically run on). So, you could just skip that test entirely and always just do the mxSetPr call regardless if you wanted to be a stickler about conforming code:

       mxSetPr(plhs[0], q);

Finally, I don't think q can ever return NULL from a mxRealloc call in a mex routine. I think the doc is to be interpreted as

"... a MAT standalone application or engine standalone application ..."

Best Answer

Related Solutions

MATLAB: Can A = A + B’*B be sped up somehow? It is seriously bottlenecking the for-loop

MATLAB: How to shrink a pre-allocated array in a C-Mex

Related Question