What exactly is Batch Normalization doing? I have recently read about Batch Normalization for Deep Learning online.
Unfortunately, the notation is really inconsistent and confusing, so perhaps someone can help.
Main Question:
Let's assume we have a neural network $\mathcal{N}$ consisiting of $D_{l}$ neurons in the $l$ th hidden layer and a dataset of $N$ samples from some $d$-dimensional space, organized in a matrix $X \in \mathbb{R}^{N \times d}$.
Then, the outputs (activations) of the $(l-1)$th layer are given by a matrix $H_{l-1} \in \mathbb{R}^{N \times D_{l-1}}$.
The input of the $i$th neuron of the next layer is hence the $i$th column of $Y_{l} = H_{l-1}W_{l} + \theta_{l}$ with each entry corresponding to one instance in the dataset.
Now what exactly is being normalized?
I would assume the following:
$$\hat Y_{l}^{ij} = \gamma \cdot \frac{Y_{l}^{ij} - \mu_{j}}{\sigma_{j}} + \beta$$
for $\mu_{j} = \frac{1}{N} \cdot \sum_{i=1}^{N} Y_{l}^{ij}$ and $\sigma_{j}$ accordingly.
Is this correct?
Finally, do the scale and offset parameters $\gamma$ and $\beta$ depend on $j$ also, or are they computed for each neuron individually?
Please can someone just give me a formula...
Bonus:
If someone can explain how this arithmetic is extended to the case if our input is a tensor used in image classification where $\dim(X) = (N,C,W,L)$ where $C$ is the number of channels, I would be very grateful, but if not I am also happy.
I usually post on the mathematics-stackexchange but this really seemed to be more appropriate here.
 A: I'm guessing that by $j$ you mean the index of the batch, i.e. $j=1$ means 1st batch, right?
What is happening is that each column $i$ gets normalized to zero mean and unit standard deviation and then shifted and scaled by $\beta$ and $\gamma$, accordingly.
This means that since you have $D_{l-1}$ columns in $H_{l-1}$:
$\mu, \sigma, \beta$ and $\gamma$ all will be vectors with $D_{l-1}$ dimensions, the latter two of which are trainable.
Thus the batch normalization operation with input $Y_{l}^{ij}$ and output $\hat Y_{l}^{ij}$ would look like this.
$$\hat Y_{l}^{ij} = \gamma_j \cdot \frac{Y_{l}^{ij} - \mu_{j}}{\sigma_{j}} + \beta_j$$
In image datasets where you have a shape of $(N, H, W, C)$, where $C$ is the number of channels, each of the variables of barchnorm  $\mu, \sigma, \beta$ and $\gamma$ would have $C$ dimensions.

We can user keras to confirm this on our own.
1)  Tabular data
import tensorflow as tf  # requires tensorflow >= 2.0.0

inp = tf.keras.layers.Input((30,))  # 30 columns (irrelevant to BN)
x = tf.keras.layers.Dense(50)(inp)  # 50 neurons on the first hidden layer
bn = tf.keras.layers.BatchNormalization()(x)  # add batchnorm after hidden layer
out = tf.keras.layers.Dense(5)(bn)  # 5 classes (irrelevant to BN)

model = tf.keras.models.Model(inp, out)
model.summary()

This will print the following:
Layer (type)                 Output Shape              Param #   
=================================================================
input_3 (InputLayer)         [(None, 30)]              0         
_________________________________________________________________
dense_2 (Dense)              (None, 50)                1550      
_________________________________________________________________
batch_normalization_2 (Batch (None, 50)                200       
_________________________________________________________________
dense_4 (Dense)              (None, 5)                 255       
=================================================================
Total params: 2,005
Trainable params: 1,905
Non-trainable params: 100
_________________________________________________________________

What interests us is the $200$ parameters that batchnorm has. Why $200$? Because there are $4$ variables (i.e. $\mu, \sigma, \beta$ and $\gamma$), each having $50$ dimensions (i.e. as many as the neurons of the previous layer).
2)  Image data
Let's do the same thing on a CNN for image classification.
inp = tf.keras.layers.Input((100, 200, 3))  # height=100px, width=200px, channels=3
c = tf.keras.layers.Conv2D(30, (4, 4), padding='same')(inp)  # same padding to keep the same height/width
bn = tf.keras.layers.BatchNormalization()(c)  # add batchnorm after conv
fl = tf.keras.layers.Flatten()(bn) 
out = tf.keras.layers.Dense(10)(fl)  # 10 classes

model = tf.keras.models.Model(inp, out)
model.summary()

This will print the following:
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
input_1 (InputLayer)         [(None, 100, 200, 3)]     0         
_________________________________________________________________
conv2d (Conv2D)              (None, 100, 200, 30)      1470      
_________________________________________________________________
batch_normalization (BatchNo (None, 100, 200, 30)      120       
_________________________________________________________________
flatten (Flatten)            (None, 600000)            0         
_________________________________________________________________
dense (Dense)                (None, 10)                6000010   
=================================================================
Total params: 6,001,600
Trainable params: 6,001,540
Non-trainable params: 60
_________________________________________________________________

Again we are interested in the $120$ parameters of batchnorm. Why $120$? Because each of the $4$ variables has $C=30$ dimensions.
