Hello KFrank,
You are absolutely correct. My optimizer was Adam with weight decay, and using SGD without weight decay gives the “expected” behaviour:
#optim = Adam(params, lr, weight_decay=1e-6)
optim = SGD(params, lr)
...

In[2]: optim.param_groups[0]["params"][0]
Out[2]: 
Parameter containing:
tensor([[ 0.0913,  0.0789,  0.0728,  ..., -0.0352,  0.0444, -0.0342],
        [ 0.0287,  0.0179,  0.0200,  ..., -0.0781, -0.0859, -0.0117],
        [-0.0233,  0.0586, -0.0877,  ..., -0.0872,  0.0094, -0.0664],
        ...,
        [-0.0286,  0.0854, -0.0020,  ..., -0.0824,  0.0082,  0.0802],
        [-0.0504, -0.0254,  0.0182,  ..., -0.0826,  0.0746,  0.0196],
        [-0.0585, -0.0425, -0.0545,  ...,  0.0924,  0.0455, -0.0788]],
       device='cuda:0', requires_grad=True)
In[3]: optim.zero_grad()
In[4]: optim.param_groups[0]["params"][0].grad
Out[4]: 
tensor([[0., 0., 0.,  ..., 0., 0., 0.],
        [0., 0., 0.,  ..., 0., 0., 0.],
        [0., 0., 0.,  ..., 0., 0., 0.],
        ...,
        [0., 0., 0.,  ..., 0., 0., 0.],
        [0., 0., 0.,  ..., 0., 0., 0.],
        [0., 0., 0.,  ..., 0., 0., 0.]], device='cuda:0')
In[5]: optim.step()
In[6]: optim.param_groups[0]["params"][0]
Out[6]: 
Parameter containing:
tensor([[ 0.0913,  0.0789,  0.0728,  ..., -0.0352,  0.0444, -0.0342],
        [ 0.0287,  0.0179,  0.0200,  ..., -0.0781, -0.0859, -0.0117],
        [-0.0233,  0.0586, -0.0877,  ..., -0.0872,  0.0094, -0.0664],
        ...,
        [-0.0286,  0.0854, -0.0020,  ..., -0.0824,  0.0082,  0.0802],
        [-0.0504, -0.0254,  0.0182,  ..., -0.0826,  0.0746,  0.0196],
        [-0.0585, -0.0425, -0.0545,  ...,  0.0924,  0.0455, -0.0788]],
       device='cuda:0', requires_grad=True)
In[7]: optim.param_groups[0]["params"][0].grad
Out[7]: 
tensor([[0., 0., 0.,  ..., 0., 0., 0.],
        [0., 0., 0.,  ..., 0., 0., 0.],
        [0., 0., 0.,  ..., 0., 0., 0.],
        ...,
        [0., 0., 0.,  ..., 0., 0., 0.],
        [0., 0., 0.,  ..., 0., 0., 0.],
        [0., 0., 0.,  ..., 0., 0., 0.]], device='cuda:0')

Thank you!