Ускорение нейронной сети GPU

Question

Я пытаюсь реализовать нейронную сеть для работы на GPU с использованием библиотек Thrust и CUBLAS, но у меня много проблем с тем, чтобы она работала быстрее, чем наша текущая многопоточная и векторизованная реализация ЦП. В сети есть один скрытый слой с логистическими единицами и выходной слой с линейными единицами, и вот код для этого:

// Functor to add bias before computing logistic
template <typename T>
struct bias_logistic_f {
        __host__ __device__
        T operator()(const T& x, const T& y) const {
                return 1/(1+exp(-(x+y)));
        }
};
bias_logistic_f bias_logistic();

// Thrust vectors for input/hidden/output units
thrust::device_vector<FLT> batch(batch_rows*ndim);
thrust::device_vector<FLT> hid(batch_rows*nhid);
thrust::device_vector<FLT> gpu_code(ndata*ncode);

// ...Load data and network weights...

// Multiply input (batch) by weights (vis2hid)
// Our matrices are stored row-major, but BLAS wants column-major,
// so pretend they're transposed and compute hid' = vis2hid' * batch'
cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, nhid, batch_rows, ndim,
            &alpha, thrust::raw_pointer_cast(&vis2hid[0]), nhid,
                    thrust::raw_pointer_cast(&batch[0]), ndim,
             &beta, thrust::raw_pointer_cast(&hid[0]), nhid);

// Add hidbiases to hid and compute logistic
thrust::transform(hid.begin(), hid.end(), hidbiases.begin(), hid.begin(),
                  bias_logistic);

// Multiply hid by weights (hid2code)
cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, ncode, batch_rows, nhid,
            &alpha, thrust::raw_pointer_cast(&hid2code[0]), ncode,
                    thrust::raw_pointer_cast(&hid[0]), nhid,
             &beta, thrust::raw_pointer_cast(&gpu_code[b*batch_rows*ncode]), ncode);

// Add codebiases
thrust::transform(gpu_code.begin() + b*batch_rows*ncode, gpu_code.begin() + (b+1)*batch_rows*ncode,
                  codebiases.begin(), gpu_code.begin() + b*batch_rows*ncode,
                  thrust::plus<FLT>());

Наши входные данные представляют собой разреженную матрицу с около 150000 строк и 6500 столбцов, в среднем около 100 ненулевых элементов на строку. Это слишком велико, чтобы хранить полную матрицу в виде плотной матрицы на графическом процессоре, поэтому я делаю цикл по разреженной матрице, расширяя партии по 1000 строк каждая для ввода в нейронную сеть:

for(int b=0; b<nbatch; ++b) {
    // Zero out batch b
    thrust::fill(batch.begin(), batch.end(), 0.0f);
    // batch_val contains the non-zero values for the current batch, batch_idx the indices within the batch,
    // and batch_ptr indexes into batch_val/batch_idx
    // This is like CSR format except instead of compressing rows, it's compressing submatrices of 1,000 rows
    thrust::scatter(batch_val.begin() + batch_ptr[b],
                    batch_val.begin() + batch_ptr[b+1],
                    batch_idx.begin() + batch_ptr[b],
                    batch.begin());

    // ...Input batch to network (shown above)...
}

Наша реализация CPU делает то же самое, используя векторы STL. Когда я запустил оба и сравнил их время выполнения, я был удивлен, обнаружив, что код графического процессора занимает в среднем около 38 секунд для обработки наших данных, в то время как код процессора занимает всего около 27 секунд. Возможно, что это связано с тем, что графическому процессору несколько лет (Tesla C1060), а сервер - более новая 24-ядерная машина. Но все же я бы подумал, что с тысячами доступных потоков это не будет на 50% медленнее.

Любые идеи, как я могу сделать этот код работать быстрее? Я новичок в программировании на GPU, поэтому я не знаю, что я могу делать неправильно. Есть ли более эффективный способ работы с разреженными матрицами, чем то, что я здесь делаю, например, использование библиотеки CUSPARSE? Или лучше было бы вообще забыть о библиотеках высокого уровня и просто написать свои собственные ядра в CUDA, чтобы объединить шаги умножения / логистики / сложения матриц?