Here is what my understanding right now.
In official turorial 33
which we  turn on the training mode in training phase:

    for phase in ['train', 'val']:
        if phase == 'train':
            scheduler.step()
            model.train(True)  # Set model to training mode
        else:
            model.train(False)  # Set model to evaluate mode


This means even we only want to train the final layer and freeze every layer before, we still need to let BatchNorm layer get to work during training and calculate the mean/var for our own training set.
Then when the model is trained, we turn off the training mode and test it and apply the mean/var from training set to test data.
So, if we choose to calculate the fully connect layer input (or whatever layer feature to be extracted) once and save it for reuse in every epoch. The right thing to do should be turn on the training mode to calculate the mean/var for all our data (train and test). In this case, the value of the extracted features  should be equal to the tutorial approach (more easy to understand, but a little computational unefficient).
Otherwise, if we turn off the training mode to extract features, I don’t where we can get whose mean/var(s) for BatchNorm layer.