After profiling my Neural Nets' code I've realized that the method, which computes the weight changes for each arc in the network (-rate*gradient + momentum*previous_delta - decay*rate*weight), already given the gradient, is the bottleneck (55% inclusive samples).
Is there any trick to compute these values in a efficient manner?
I can't add comments; your question is incomplete as provides little information of your implementation details.
A good implementation vectorizes forward and backward propagation in the first step; then all matrix vector multiplications can be compressed into matrix matrix multiplication by training with mini batches.
If you are done with those steps and your gradients passed with less than 10^-4 tolerance you should check out ?gemm_ [matrix * matrix] calls in fast matrix algebra package which are implemented parallel and optimized to the last bit. MKL is one of those packages. Try not to implement ?gemm_ sub routines on your own, harder then it seems. Of course having large theta weights you might want to move to CUDA.
Depending on the initial criteria there may still be more room for optimization; providing the data set and weights are large enough you can go parallel; or even explore distributed weight updates.
good luck
