I cannot reproduce the increased memory usage with my manual implementation and comparing it via:
bn = MyBatchNorm2d(3, affine=True).cuda()
# bn = nn.BatchNorm2d(3, affine=True).cuda() # switch to this in another run
for _ in range(10):
    x = torch.randn(16, 3, 224, 224, device='cuda')
    out = bn(x)
    out.mean().backward()
    print(torch.cuda.memory_allocated()/1024**2)
    bn.zero_grad()

and both print a memory allocation of ~18MB.