I am confused on the application/calculation of BN statistics in hidden layers in deep nets. In general, any input can be normalized and the mean and variance statistics across the mini-batch can be computed at the first layer. This scaled input is then fed to the network. However, the outputs of deeper layer are in general of different dimensions than the input. How are the activations in deeper layers scaled? How are the means and variance computed for a random deeper layer in the network? Its not clear to me how the changing matrix dimensions is taken into account.
Batch normalization is not a method of standardizing the inputs; rather, it standardizes the activations of the hidden units. For example, if your network has $l$ layers, batch normalization could be applied separately to each layer by subtracting the mean from and dividing by the standard deviation of the activations of the hidden units in each layer (during forward propagation).
Moreover, in practice, additional terms are used to shift the mean and scale the standard deviation of the units in a layer to arbitrary values, which are learned through backprop. These values may be different from those that result from the outputs of the hidden units in a layer because they are learned directly rather than implicitly as a function of those hidden units' activations.
Broadly speaking, batch normalization is important because the gradients of the units can be very small or very large, depending on the nonlinearity the network is modeling. This makes it difficult to choose a common learning rate for gradient descent that can be applied to all the units because layers may affect each other exponentially. Normalization is useful because it scales the units outputs to a certain scale which makes gradient descent faster and, in a sense, more efficient.
For more information, please see Deep Learning by Goodfellow et al. (2016) or Ioffe and Szegedy (2015, referenced therein).
Having different dimensions at different layers doesn't have anything to do with normalization or the mean and sigma of either weights or activation functions at each layer. Activation functions at each node of a hidden layer are typically either the logistic() or tanh() function, whose outputs are constrained to boundary conditions. You are mean-zero standardizing or normalizing the input features to range [0,1] so that they are constrained as well.
Exactly how the mean and sigma are determined is highly specific to a network model, and there is an infinitely large number of ways to define means and variance in an ANN, since there are so many parameters. The main point is that learning via an ANN almost always needs to be constrained, or the solution won't converge. If the values of the updated connection weights between nodes start wildly shooting out of their historical range over the history of iterations, then the model will depart and not converge. In this case, error won't decrease it will blow up. This can happen when the learning rate is off, the momentum is too large, or weight decay is off.
You might lookup parameter "tuning" via use of a grid search before your models starts learning. This merely involves changing starting values and looping through your data several times, to see which values resulted in the greatest reduction in error. Then start your full training using the tuned parameters.