Q:

Speeding up Matlab optimizations by CUDA

I was wondering if I can improve the speed of my simulations with CUDA in Matlab by simply breaking up the loop into GPU-threads.
So far I have done the loops with mclapply in R, but I was wondering, if I could speed it up in Matlab with CUDA? (Nvidia Gtx 650 ti boost)
Here is a simplified example of what I'm doing:
I have a dataset of dim 2000x9 and I want to run at least 250 loops with moving data window (ca. 1500 rows), these loops are independent thus perfect for parallel computation.
Example of loop: Take the first 1500 rows. Some magic with the whole dataset (1500x9). Calculate univariate function on each column (do 9times operation on 1500x1) and then run minimization problem with a certain loss function for each column (do 9times operation on 1500x1).
(inbetween I'm also working with all columns, so I cannot separate the columns in different threads)
My idea:
If I would be running each loop in different GPU-thread? Does it make sense to do such a hard computational problem in single GPU core? (for illustration, 1 loop takes ca. 300seconds on i7 3770k on 1 core)
I know that single GPU threads are much slower compared to CPU, but if all loops will be running simultaneously? Especially since each loop takes so long once its given data, so feeding the data to GPU will take only insignificant portion of overall time.
Important note: I'm a very lousy "programmer", and what I plan is just splitting the for loop (i in 1:250) in function of variable i and passing this to GPU.

A:

I see three main issues related to the fact that each thread should perform a whole minimization process (which kind of Matlab minimization routine are you currently using? fminunc? fminsearch? minFunc?). 

Minimizations can be demanding in terms of needed temporary variables. This could limit the performance of your algorithm due to the need of storing temporaries and transacting them, for example, with global memory, depending of course on how do you mean to implement it. 
You should also carefully take into account thread synchronization, since the time needed for accomplishing the minimization process could change from thread to thread.
Matlab has very effective optimization routines, whose performance are generally difficult (but, of course, not impossible) to replicate by a custom implementation. In my experience, Matlab's fminunc is more effective than the Broyden-Fletcher-Goldfarb-Shanno equivalent routine provided by the NAG. So, if you are trying to translate one of the above optimization routines, then you could end up by a less satisfactory result.

I have faced many optimization problems using Matlab accelerated with CUDA and my "gold rule" is to use one of the Matlab's optimization routines and accelerate the solution of the direct problem (calculation of the functional) and the functional gradient by purposely written CUDA codes interfaced with Matlab by mex-files. Take into account that especially the gradient needs to (and can) be accelerated since the calculation of the functional derivative by finite differences are independent and require the invokation of as many functional computation routines as is the number of optimization parameters.
EDIT
Suppose that I have to optimize the objective functional objfun. What I'm doing is to code objfun in CUDA with a mex-file interface, compiling it by nvcc and then linking it under Matlab.
As I'm using Matlab 2010, the CUDA function is compiled by nvcc and transformed in C++ code by the command
system(sprintf('nvcc -I"%s/extern/include" --cuda "mexfun.cu" --output-file "mexfun.cpp"', matlabroot));

and then linked to Matlab by
mex -I/opt/cuda/include -L/opt/cuda/lib -lcudart mexfun.cpp

as suggested in Compiling CUDA C/C++ mex code under linux.
Then, when using, for example, fminunc(@mexfun,...), Matlab will optimize the objective functional and each evaluation of it will be performed (and thus accelerated) on the GPU. When analytically available, I'm also coding the gradient calculation by the same approach, since finite differences used to evaluate the gradient can significantly slow down the whole optimization process.
For Matlab 2013 and Windows systems, see Creating mex files from CUDA code.
EDIT Structure of mexfun.cu (objective function)
// Do not change the function name (`mexFunction`) and the function arguments (`nlhs`, `plhs`, ...). 
void mexFunction(int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[])

{
    /* Maps Matlab's pointers to the input variables to CUDA pointers */
    double* input_1     = mxGetPr(prhs[0]);
    double* input_2     = mxGetPr(prhs[1]);

    /* Recovers the size of the input matrices */
    int dimx = mxGetN(prhs[0]);
    ...         
    int dimu = mxGetM(prhs[3]);         

    /* Memory allocations on the host */
    cuDoubleComplex* hfoo = (cuDoubleComplex *)malloc(sizeof(cuDoubleComplex)*dimx);
    ...

   /* Memory allocations on the device */
   cuDoubleComplex* dfoo; cudaMalloc((void*)&d_Kernel_Matrix,dimx*sizeof(cuDoubleComplex));
   ...

  /* Memory transfer from host to device */
  cudaMemcpy(dfoo,hfoo,dimx*sizeof(cuDoubleComplex),cudaMemcpyHostToDevice);
  ....

  /* Kernel launch */
  dim3 dimBlock(BLOCK_SIZE_X,BLOCK_SIZE_Y);
  Kernel_To_Be_Launched <<<dimGrid,dimBlock >>>(hfoo,dfoo,dimx);

 /* Copy the results from device to host */ cudaMemcpy(hfoo,dfoo,dimx*sizeof(cuDoubleComplex),cudaMemcpyDeviceToHost);

 /* Passing the output matrices to MATLAB */
 plhs[0] = mxCreateDoubleMatrix(1,dimu,mxCOMPLEX);
 double* hfoo_re = mxGetPr(plhs[0]);
 double* hfoo_im = mxGetPi(plhs[0]);

 /* Freeing host memory */
 free(hfoo);
 ...

 /* Freeing device memory */
 cudaFree(dfoo);

}

