Ah, I see now. The derivative for LogSoftMax was transposed from how I was thinking about it conceptually:
dlogits_dnew_relu_feat = torch.Tensor([[[1-sm[0][0], -sm[0][0], -sm[0][0]], [-sm[0][1], 1-sm[0][1], -sm[0][1]], [-sm[0][2], -sm[0][2], 1-sm[0][2]]]])

should  be
dlogits_dnew_relu_feat = torch.Tensor([[[1-sm[0][0], -sm[0][1], -sm[0][2]], [-sm[0][0], 1-sm[0][1], -sm[0][2]], [-sm[0][0], -sm[0][1], 1-sm[0][2]]]])

and now you can do dloss_dlogits @ dlogits_dnew_relu_feat with no need for squeeze/unsqueeze