For the simple repro example, the Python equivalent gives the same output as BatchNorm2d in libtorch 1.4:
>>> import torch
>>> torch.manual_seed(1)
<torch._C.Generator object at 0x7f3d14d2ddb0>
>>> input = torch.rand(1, 3, 4, 4)
>>> torch.manual_seed(1)
<torch._C.Generator object at 0x7f3d14d2ddb0>
>>> torch.nn.BatchNorm2d(3)(input)
tensor([[[[ 1.1911, -1.1530, -0.5465,  1.0787],
          [-2.3783,  1.3981, -0.5756,  1.1752],
          [ 0.2692, -0.3715,  0.6082,  0.0494],
          [ 0.8235, -1.0264, -0.2501, -0.2921]],
 
         [[ 0.0975, -0.1822,  1.4672,  0.4098],
          [-0.8742, -1.3091, -0.4895, -0.9849],
          [-0.7766, -0.0848,  0.9450,  0.8462],
          [-2.0108,  0.9897,  0.3500,  1.6069]],
 
         [[ 1.0529, -1.4680, -1.5317, -0.7800],
          [ 1.4028, -0.2734,  0.6805, -0.7517],
          [ 1.5682, -0.6850,  1.1971,  0.0129],
          [-0.8514,  0.8186, -0.8578,  0.4658]]]],
       grad_fn=<NativeBatchNormBackward>)

For model convergence in general, have you tried training the equivalent Python model and checking that it converges?