I have checked the result and found this phenomenon:
if I use “version 1” instead of  MyBNFunc, everything works well and it trained successfully.
but if I turned to use the current version, the training failed but the forward output matched exactly.
however, I implemented my function according to the BatchNorm paper, and cannot find what is wrong with it.
class MyBNFunc(Function):
    @staticmethod
    def forward(ctx, input, avg, var, gamma, beta, eps):
        B, C, H, W = input.shape
        avg = avg[None, :, None, None]
        var = var[None, :, None, None]
        ctx.avg = avg
        ctx.var = var
        ctx.eps = eps
        ctx.B = B
        output = input - avg
        output = output / (torch.sqrt(var + eps))
        ctx.save_for_backward(input, gamma, beta, output)
        scaled_output = output * gamma + beta
        return scaled_output
    
    @staticmethod
    def backward(ctx, grad_output):
        input, gamma, beta, output = ctx.saved_tensors
        avg = ctx.avg
        var = ctx.var
        eps = ctx.eps
        B = ctx.B
        dL_dxi_hat = grad_output * gamma
        dL_dvar = (-0.5 * dL_dxi_hat * (input - avg) / ((var + eps) ** 1.5)).sum((0, 2, 3), keepdim=True)
        dL_davg = (-1.0 / torch.sqrt(var + eps) * dL_dxi_hat).sum((0, 2, 3), keepdim=True) + dL_dvar * (-2.0 * (input - avg)).sum((0, 2, 3), keepdim=True) / B
        dL_dxi = dL_dxi_hat / torch.sqrt(var + eps) + 2.0 * dL_dvar * (input - avg) / B + dL_davg / B # dL_dxi_hat / sqrt()
        dL_dgamma = (grad_output * output).sum((0, 2, 3), keepdim=True) #(grad_output * output).sum((0, 2, 3), keepdim=True)
        dL_dbeta = (grad_output).sum((0, 2, 3), keepdim=True)
        return dL_dxi, None, None, dL_dgamma, dL_dbeta, None
    
    
class MyBN1(nn.BatchNorm2d):
    def __init__(self, num_features, eps=1e-5, momentum=0.1,
                 affine=True, track_running_stats=True):
        super(MyBN1, self).__init__(
            num_features, eps, momentum, affine, track_running_stats)

    def forward(self, input):
        self._check_input_dim(input)

        exponential_average_factor = 0.0

        if self.training and self.track_running_stats:
            if self.num_batches_tracked is not None:
                self.num_batches_tracked += 1
                if self.momentum is None:  # use cumulative moving average
                    exponential_average_factor = 1.0 / float(self.num_batches_tracked)
                else:  # use exponential moving average
                    exponential_average_factor = self.momentum

        # calculate running estimates
        if self.training:
            mean = input.mean([0, 2, 3])
            # use biased var in train
            var = input.var([0, 2, 3], unbiased=False)
            n = input.numel() / input.size(1)
            with torch.no_grad():
                self.running_mean = exponential_average_factor * mean\
                    + (1 - exponential_average_factor) * self.running_mean
                # update running_var with unbiased var
                self.running_var = exponential_average_factor * var * n / (n - 1)\
                    + (1 - exponential_average_factor) * self.running_var
        else:
            mean = self.running_mean
            var = self.running_var
        weight = self.weight[None, :, None, None]
        bias = self.bias[None, :, None, None]

        output = MyBNFunc.apply(input, mean, var, weight, bias, self.eps)
        '''
        VERSION 1
        input = (input - mean[None, :, None, None]) / (torch.sqrt(var[None, :, None, None] + self.eps))
        if self.affine:
            input = input * self.weight[None, :, None, None] + self.bias[None, :, None, None]
        '''

        return output