Note that R[k, j] = R[k, j] - regularization * R[k, j] is equivalent to R[k, j] = R[k, j] * (1 - regularization). Instead of going through all elements of R[:, j] and doing this operation, you can store the scale of R[:, j] in R_scales[j]. This way the cost of the regularization is just one multiplication. Then when you compute the estimatation, you do your dot product and multiply the result by R_scales[j]. Of course the same goes for L. 

Grep for w_scale in sgd_fast.pyx. This can speed up the code by 50% in my experience.