Hi,
Weight decay in the optimizers is computed “in closed form” by simply doing
grad += weight_decay * weight

(done here for sgd for example).
This will be more numerically stable and will have less floating point errors because it avoid doing the whole backward pass.
So the two are expected to be different but they should not lead to big change in convergence unless your net is big enough that the sum of the squared weights is big enough to have floating point precision problem.
Given what you say, this might actually be happening.
You can quickly check that by trying to do the accumulation in double:
wd_double = 0.
wd_single = 0.
for param in net.parameters():
  wd_double += weight_decay * 0.5 * torch.pow(param.double(), 2).sum()
  wd_single += weight_decay * 0.5 * torch.pow(param, 2).sum()
print("Double precision: ", wd_double)
print("Single precision: ", wd_single)
print("Difference: ", wd_double - wd_single)

I would actually be interested to see what your net is and what this returns?