            for (int x = xi_l; x <= xi_r; ++x) {
              if (!between(x, 0, width-1)) continue;
              int address = inputImages_strideBatch * b + inputImages_strideHeight * y + inputImages_strideWidth * x;
-             grad_yf += 2 * gradOutValue_fg * inputImages_data[address + t]
+             grad_yf += gradOutValue_fg * inputImages_data[address + t]
                         * weights[id_point]
                         * (bias_y[id_point] - bias_y_weighted / weight_sum)
-                        / weight_sum;
-             grad_xf += 2 * gradOutValue_fg * inputImages_data[address + t]
+                        / weight_sum / 2;
+             grad_xf += gradOutValue_fg * inputImages_data[address + t]
                         * weights[id_point]
                         * (bias_x[id_point] - bias_x_weighted / weight_sum)
-                        / weight_sum;
+                        / weight_sum / 2;
              ++id_point;
            }
          }

           if (!between(x, 0, width-1)) continue;
           bias_x[id_point] = (x - xf_s); // remember multiply 2 below since it is derivative of square
           bias_y[id_point] = (y - yf_s); // remember multiply 2 below since it is derivative of square
           float weight = expf(-bias_x[id_point] * bias_x[id_point] - bias_y[id_point] * bias_y[id_point]);
           bias_x_weighted += weight * bias_x[id_point];
           bias_y_weighted += weight * bias_y[id_point];
           weight_sum += weight;

static int cunn_SubSamplerBHWD_updateOutput(lua_State *L)
{
  THCState *state = getCutorchState(L);

  THCudaTensor *inputImages = (THCudaTensor *)luaT_checkudata(L, 2, "torch.CudaTensor");
  THCudaTensor *grids = (THCudaTensor *)luaT_checkudata(L, 3, "torch.CudaTensor");
  THCudaTensor *masks = (THCudaTensor *)luaT_checkudata(L, 4, "torch.CudaTensor");
