Introduction

Here is a high-level summary of the material so far. Using linear classification as an example, we introduced softmax and SVM classifiers and discussed neural networks as a kind of score function used to train them. We examined the structure of a perceptron, which models a biological neuron, and the meaning of each operation. We then moved beyond simple logical structures to deep neural networks with multiple nonlinear layers and saw how chain-rule optimization makes them universal approximators for real-world tasks involving complex functions. Greater depth increases representational power—the ability to express complex functions—but also causes overfitting, such as growing weight variance or excessive adaptation to the training set. We therefore introduced regularization methods including dropout and $L_2$ regularization.
This post turns from the architecture and significance of deep networks to practical training concerns: data preprocessing, weight initialization, and the benefits of regularizers such as batch normalization.

Data Preprocessing

Let us apply several transformations to data $X$. Assume $X$ is a batch matrix $X \in \mathbb{R}^{N \times D}$, where $N$ is the number of samples and $D$ their dimensionality. For ten RGB images of size $3 \times 32 \times 32$, for example, $N=10$ and $D=3072$.

Original data

First consider the original data, meaning the data before any preprocessing, and assume it has the following distribution.

For simplicity, assume two-dimensional data, with each point in the coordinate plane representing a sample.

Mean subtraction

One preprocessing method is mean subtraction, which removes bias in the data distribution by subtracting the mean of each dimension.

original_data = X #(assume that X is numpy array with shape N*D)
zero_centered_data = X - np.mean(X, axis = 0)

For readers unfamiliar with numpy, np.mean(_, axis=0) computes the mean along the array's zeroth axis. For two-dimensional data, this means subtracting the mean of each column vector, each in $\mathbb{R}^{N \times 1}$. Removing the bias with each column's mean moves the center of the original distribution to the origin $(0,0)$.

Normalization

Each dimension can be understood as a feature. If different dimensions have different variances, as above, their relative importance or effective learning rates during optimization may differ. Preprocessing should make their scales more similar.

original_data = X #(assume that X is numpy array with shape N*D)
zero_centered_data = X - np.mean(X, axis = 0)
normalized_data = zero_centered_data / np.std(X, axis = 0)

PCA and Whitening

The preceding methods normalize dataset samples. Feature-vector columns generally remain correlated, as in $y=ax$. PCA and whitening remove this correlation.

original_data = X #(assume that X is numpy array with shape N*D)
zero_centered_data = X - np.mean(X, axis = 0)
covariance_matrix = zero_centered_data.T.dot(zero_centered_data)/X.shape[0]

The process is as follows. Center data $X \in \mathbb{R}^{N \times D}$ at zero, take the centered matrix's inner product with itself, and divide by sample count $N$ to obtain the covariance matrix.

$$Cov(X) = \frac{\left(X-\mu(X) \right)^\top \left(X-\mu(X) \right)}{N}$$

Each covariance-matrix element represents the relationship between features $i$ and $j$; the example yields a $2 \times 2$ matrix. A covariance matrix is automatically symmetric. Its diagonal contains each feature's autocorrelation, which mathematically equals its variance. Singular value decomposition (SVD) of the covariance matrix extracts eigenvectors $U$, eigenvalues $V$, and singular values $S$.

U, S, V = np.linalg.svd(covariance_matrix)

SVD is useful because covariance-matrix elements express correlations between features. Orthonormal matrix $U$ serves as a projection basis that removes correlations from the mean-centered original data.

X_dr = zero_centered_data.dot(U)

The dot product rotates the axes according to orthonormal basis $U$. Since the original feature vectors were correlated, rotating the axes removes those correlations—a decorrelation step.

PCA stands for Principal Component Analysis. Rather than using every feature basis extracted by SVD, it retains only the important ones. If eigenbases are sorted by descending eigenvalue, we might keep the 100 most important features among $D$. PCA is unnecessary for the two-dimensional example, but selecting meaningful features becomes important as dimensionality grows. The curse of dimensionality arises because the data manifold often lies not throughout the ambient space, but on a particular surface such as a hyperplane within it.

Xrot_reduced = np.dot(X, U[:,:100]) # Xrot_reduced becomes [N x 100]

The Swiss roll dataset is three-dimensional, but its distribution forms a rolled-up plane. Projecting it onto the plane that best represents the data distribution, rather than using the raw 3D space, is an example of PCA.
Reducing the earlier example from $N \times D$ to $N \times 100$ gives fewer, decorrelated features. After projecting the data onto the eigenbasis, dividing each dimension by its eigenvalue normalizes by the variance of samples along that basis. This is whitening. $X$ has already been rotated into the basis, and the eigenvalues contain the variance of each feature because the corresponding $\sigma$ values determine the diagonal elements of the covariance matrix.

# whiten the data:
# divide by the eigenvalues (which are square roots of the singular values)
Xwhite = Xrot / np.sqrt(S + 1e-5)

The name “whitening” comes from producing a white-noise-like distribution with zero mean and an identity covariance matrix.

Weight initialization

Unlike instance-based learning such as K-nearest neighbors, in which training instances directly participate in test predictions, model-based learning learns by updating weight parameters determined by the model architecture. Performance therefore depends on how weights are initialized. What is the most effective method?

Zero initialization

The first idea is to initialize every weight to $0$. Although we cannot predict every optimized parameter value, the law of large numbers suggests that in a deep network with many parameters, roughly half the weights will become positive and half negative. It may therefore seem reasonable to start at mean zero and let some weights move positive and others negative, but this is wrong. If every neuron emits the same output, backpropagation updates them identically. Initializing every weight to the same value leaves the weight matrix symmetric. The following perceptron code provides a simple example.

import numpy as np

# sigmoid as an activation function
def sigmoid(x):
    return 1 / (1 +np.exp(-x))

# input : 2 * 3 [N * D]
input = np.array([[-1, 0, 1], [-1, 1, 0]])

# target : 2 * 2 [N * out]
target = np.array([[0, 1], [1, 0]])

# Weight 1 [3 * 4]
W1 = np.zeros((3, 4))

# Weight 2 [4 * 2]
W2 = np.zeros((4, 2))

# feed forward
hidden = sigmoid(input.dot(W1))
out = hidden.dot(W2)

# calculate difference
diff = target - out

# backpropagation
W2 -= hidden.T.dot(diff)
W1 -= input.T.dot((hidden*(1-hidden)))

# print results
print(W1, W2)

Because the code is hard to read, the same computation can be written as follows.

$$\begin{aligned} \text{input} =& \begin{bmatrix} -1 & 0 & 1 \newline -1 & 1 & 0 \end{bmatrix},~\text{target} = \begin{bmatrix} 0 & 1 \newline 1 & 0 \end{bmatrix} \newline \text{output} =& \sigma \left( input \odot W1 \right) \odot W2 \newline \text{diff} =& \text{target} - \text{output} \end{aligned}$$

Inspecting the weight parameters after running the code gives

$$W_1 = \begin{bmatrix} 0.5 & 0.5 & 0.5 & 0.5 \newline -0.25 & -0.25 & -0.25 & -0.25 \newline -0.25 & -0.25 & -0.25 & -0.25 \end{bmatrix},~W_2 = \begin{bmatrix} -0.5 & -0.5 \newline -0.5 & -0.5 \newline -0.5 & -0.5 \newline -0.5 & -0.5 \end{bmatrix}$$

Every value in each row is updated identically. The number of weight parameters effectively determines the number of usable nodes in the neural network, but such identical updates make the extra parameters pointless. In other words, even with many parameters, representational power collapses. This problem occurs whenever all parameters share one initial value, not only when that value is zero.

import numpy as np

# sigmoid as an activation function
def sigmoid(x):
    return 1 / (1 +np.exp(-x))

input = np.array([[-1, 0, 1], [-1, 1, 0]])
target = np.array([[0, 1], [1, 0]])

value1 = 0.01
value2 = 0.09

W1 = value1 * np.ones((3, 4))
W2 = value2 * np.ones((4, 2))

# feed forward
hidden = sigmoid(input.dot(W1))
out = hidden.dot(W2)

# calculate difference
diff = target - out

# backpropagation
W2 -= hidden.T.dot(diff)
W1 -= input.T.dot((hidden*(1-hidden)))

# print results
print(W1, W2)

Here, W1 and W2 were each initialized to identical values before training. The result is as follows.

$$W_1 = \begin{bmatrix} 0.51 & 0.51 & 0.51 & 0.51 \newline -0.24 & -0.24 & -0.24 & -0.24 \newline -0.24 & -0.24 & -0.24 & -0.24 \end{bmatrix},~W_2 = \begin{bmatrix} -0.23 & -0.23 \newline -0.23 & -0.23 \newline -0.23 & -0.23 \newline -0.23 & -0.23 \end{bmatrix}$$

Small random numbers

To avoid overfitting and the problem above, weights cannot be initialized to the same value. The next option is random initialization.

W = 0.01* np.random.randn(D, H)

Here, $H$ is the hidden layer's output dimension. randn generates a $D \times H$ matrix of normal random values, so the code produces weight matrix $W \in \mathbb{R}^{D \times H}$ from a Gaussian with standard deviation $0.01$.

Calibrating the variances with $\frac{1}{\sqrt{N}}$

The problem with this method is that output variance can grow with the number of inputs. Instead of choosing an arbitrary Gaussian standard deviation such as 0.01, we want each neuron's output distribution to have variance—or standard deviation—$1$. If a layer has $n$ inputs, initialize its neuron weights as

w = np.random.randn(n) / sqrt(n)

Doing so gives every neuron in the network the same output distribution and can increase convergence speed.

$$\begin{aligned} \text{Var}(s) =& \text{Var} \left( \sum_i^n w_i x_i \right) \newline =& \sum_i^n \text{Var} (w_i x_i) \newline =& \sum_i^n \left( E(w_i) \right)^2 \text{Var} (x_i) + \left( E(x_i) \right)^2 \text{Var}(w_i) + \text{Var} (x_i) \text{Var}(w_i) \newline =& \sum_i^n \text{Var} (x_i) \text{Var}(w_i) \newline =& (n \text{Var}(w)) \text{Var} (x) \end{aligned}$$

The derivation above proves this concept.

Sparse initialization

Solving the calibration problem—different output distributions across neurons—might suggest zero-initializing every weight matrix, but then all weights learn identically. Instead, fix the number of active neurons and compute only weights sampled from random Gaussian noise—the small, uncalibrated random values seen earlier. Fixing the amount of computation keeps the output distribution consistent.

Initializing bias

Biases may safely be initialized to $0$. If weight asymmetry is guaranteed, node-level asymmetry remains even when biases are identical. Since ReLU gives no gradient below zero, some practitioners initialize ReLU-neuron biases to a small positive value such as $0.01$. There is no evidence that this improves performance, however, and it can produce worse results. Zero bias initialization is therefore the most common choice.

He initialization

He et al. (reference) report that the following initialization works better for networks using ReLU. The paper contains meaningful experiments beyond weight initialization and is worth reading.

W = np.random.randn(n)*sqrt(2.0/n)

Batch normalization

Beginning with Deep Residual Learning for Image Recognition, better known as ResNet, adding batch normalization to convolutional neural networks became standard practice.

Batch normalization literally normalizes by batch. Related methods include layer, instance, and group normalization, which differ in the axes over which they normalize. Layer normalization uses the mean and standard deviation across channels within one sample; instance normalization normalizes each sample and channel; group normalization divides channels into groups within a sample and normalizes each group. This section focuses on the most common method, batch normalization.
Gradient-descent algorithms generally update parameters from mini-batches rather than individual samples. As discussed earlier, single-sample updates are noisy and cannot exploit parallel computation, making them inefficient.

Consider passing batches through a multilayer neural network. Since each batch is sampled randomly, its data distribution differs, and feature-map distributions also vary across layers. This is called internal covariate shift. It matters because it directly affects weight initialization: if batch distributions and intermediate feature-map distributions differ, the optimization cannot rely consistently on the initialization criterion. Batch normalization reduces these distributional differences by normalizing each batch, thereby reducing variability and stabilizing training.
The mean and variance of each batch are defined as follows.

$$\begin{aligned} \mu\_\text{batch} =& \frac{1}{m} \sum\_{i=1}^m x_i \newline \sigma^2\_\text{batch} =& \frac{1}{m} \left(x_i - \mu_{\text{batch}} \right)^2 \end{aligned}$$

Here, $x_i$ is sample $i$ in a batch of size $m$. The batch mean and variance are used to normalize, scale, and shift the batch.

$$\hat{x_i} = \frac{x_i - \mu\_\text{batch}}{\sqrt{\sigma\_\text{batch}^2 + \epsilon}}$$

Simple normalization alone can force a nonlinear network to remain in a linear regime, reducing representational power. Learnable parameters $\gamma$ and $\beta$ therefore scale and shift each batch-normalized pre-activation so that every layer can represent a different distribution.

$$y_i = \gamma \hat{x_i} + \beta$$

Appendix for batch normalization

Batch normalization is essential knowledge for anyone working with deep networks, yet many people apply it without fully understanding its implementation or how it works. We will therefore implement its forward and backward passes directly from the formulas. The code and explanation are based on other sources, which can be consulted directly if preferred (reference 1, reference 2).

Why and How to Use Batch Normalization

As its name implies, batch normalization normalizes input $x$ at the batch level. It commonly appears between convolutional layers and their activation functions. By preventing covariate shift—variation in feature-map and batch distributions—it enables stable training at larger learning rates and faster optimization. It should not be used indiscriminately, however. In GAN-based models, for example, batch normalization near the ends of an encoder or decoder may interfere with training. It was introduced to stabilize learning, not as a universal cure for every task. The process in the paper matches the derivation above.

$$\begin{aligned} \mu\_\text{batch} =& \frac{1}{m} \sum\_{i=1}^m x_i \newline \sigma^2\_\text{batch} =& \frac{1}{m} \left(x_i - \mu_{\text{batch}} \right)^2 \newline \hat{x_i} =& \frac{x_i - \mu\_\text{batch}}{\sqrt{\sigma\_\text{batch}^2 + \epsilon}} \newline y_i =& \gamma \hat{x_i} + \beta \end{aligned}$$

First compute the batch mean $\mu$ and variance $\sigma^2$, then normalize the samples with them. Learnable parameters $\gamma$ and $\beta$ subsequently scale and shift the result and are updated as follows.

To compute backpropagation, draw a diagram like the one above, find each local gradient by the chain rule, and multiply it by the incoming gradient. Proceeding from output to input yields gradients for every parameter along the path.

Forward propagation

N, D = x.shape

#step1: calculate mean
mu = 1./N * np.sum(x, axis = 0)

#step2: subtract mean vector of every trainings example
xmu = x - mu

#step3: following the lower branch - calculation denominator
sq = xmu ** 2

We will proceed about three steps at a time. First compute the sample mean. Since x has shape batch * dimension, averaging over axis=0 gives the batch mean. np.mean() can compute it directly instead of np.sum().

mu = np.mean(x, axis=0)

Subtract the mean from the samples. The term x - mu is reused in two computations: normalizing $x$ as below and computing the variance.

$$\begin{aligned} \hat{x_i} =& \frac{x_i - \mu\_\text{batch}}{\sqrt{\sigma\_\text{batch}^2 + \epsilon}} \newline \sigma^2\_\text{batch} =& \frac{1}{m} \left(x_i - \mu_{\text{batch}} \right)^2 \end{aligned}$$

In both expressions, the shared term is $x_i - \mu_\text{batch}$.

#step4: calculate variance
var = 1./N * np.sum(sq, axis = 0)

#step5: add eps for numerical stability, then sqrt
sqrtvar = np.sqrt(var + eps)

#step6: invert sqrtvar
ivar = 1./sqrtvar

sq is the square of x - mu, so averaging it gives the variance. As above, np.mean() can perform this operation directly.

var = np.mean(sq, axis = 0)

Adding eps in step 5 prevents division-by-zero errors. The remaining calculations produce the following output.

#step7: execute normalization
xhat = xmu * ivar

#step8: Nor the two transformation steps
gammax = gamma * xhat

#step9
out = gammax + beta

This scales and shifts xhat—the normalized x—using learnable parameters gamma and beta.

Backward propagation

$$\begin{aligned} \frac{\partial l}{\partial \hat{x_i}} =& \frac{\partial l}{\partial y_i} \cdot \gamma \newline \frac{\partial l}{\partial \sigma_B^2} =& \sum_{i=1}^m \frac{\partial l}{\partial \hat{x_i}} \cdot (x_i - \mu_B) \cdot -\frac{1}{2}(\sigma_B^2 + \epsilon)^{-3/2} \newline \frac{\partial l}{\partial \mu_B} =& \sum_{i=1}^m \frac{\partial l}{\partial \hat{x_i}} \cdot \frac{-1}{\sqrt{\sigma_B^2 + \epsilon}} \newline \frac{\partial l}{\partial x_i} =& \frac{\partial l}{\partial \hat{x_i}} \cdot \frac{1}{\sqrt{\sigma_B^2 + \epsilon}} + \frac{\partial l}{\partial \sigma_B^2} \cdot \frac{2(x_i - \mu_B)}{m} + \frac{\partial l}{\partial \mu_B} \cdot \frac{1}{m} \newline \frac{\partial l}{\partial \gamma} =& \sum_{i=1}^m \frac{\partial l}{\partial y_i} \cdot \hat{x_i} \newline \frac{\partial l}{\partial \beta} =& \sum_{i=1}^m \frac{\partial l}{\partial y_i} \end{aligned}$$

The paper summarizes backward propagation with the expression above, but implementing it is more complicated than it looks. Let us derive each gate step by step, beginning with the final gate.

$$\text{out} = \gamma \hat{x} + \beta$$

For this operation, the local gradients with respect to the addition gate's two inputs, $\gamma \hat{x}$ and $\beta$, are as follows.

dgammax = dout
dbeta = np.sum(dout, axis=0)

At first, the use of np.sum() for dbeta may seem puzzling. NumPy does not add $\beta$ at its original size; broadcasting expands it to the dimensions of gammax. Computing dbeta must reverse that broadcasting, hence the sum.

Next backpropagate through the multiplication gate. Since dgammax is the incoming gradient, the chain rule multiplies it by each local derivative to obtain the gradient flowing to each input. As discussed previously, the local derivative of a * gate uses the input from the opposite branch. The upper input is xhat and the lower is gamma, so the upper-branch gradient dxhat is dgammax * gamma, while the lower-branch gradient is dgammax * xhat.

dxhat = dgammax * gamma
dgamma = np.sum(dgammax * xhat, axis=0)

Like beta, gamma is broadcast, so np.sum() restores the proper dimensions. Backpropagation for learnable parameters gamma and beta is now complete; the remaining normalization gradients follow.

The rightmost multiplication gate is identical to the one above, so the code computes it as follows.

divar = np.sum(dxhat*xmu, axis=0)
dxmu1 = dxhat * ivar

Inverse sample variance divar was also broadcast, so it is summed with np.sum. dxmu1 equals incoming gradient dxhat multiplied by ivar. Importantly, dxmu requires two gradient paths because x - mu branched into two computations.

Backpropagation through this branch therefore computes the two gradients separately and adds them. We already have upper-path dxmu1. For the lower path, the local derivative of the inverse is $\frac{d}{dx}\left(\frac{1}{x}\right)=-\frac{1}{x^2}$, so

dsqrtvar = -1. /(sqrtvar**2) * divar

The code above follows directly. Likewise, the local gradient through the square root is $\frac{d}{dx}\sqrt{x+\epsilon}=\frac{1}{2\sqrt{x+\epsilon}}$, so

dvar = 0.5 * 1. /np.sqrt(var+eps) * dsqrtvar

This yields the expression above. The derivative of $\frac{1}{N}\sum_i x_i$ is slightly more involved because the differentiated object is a matrix, and can be defined as follows.

$$\frac{d}{dx} \left( \frac{1}{N} \sum_i x_i \right) = \frac{1}{N} \begin{pmatrix} 1 & \cdots & 1 \newline \vdots & \ddots & \vdots \newline 1 & \cdots & 1 \end{pmatrix}$$

dsq = 1. /N * np.ones((N,D)) * dvar

For the final part, $\frac{d}{dx}(x^2)=2x$, so

dxmu2 = 2 * xmu * dsq

Now that dxmu2 is available, the subtraction gate's local gradient is positive for x and negative for mu, so

dx1 = (dxmu1 + dxmu2)
dmu = -1 * np.sum(dxmu1 + dxmu2, axis=0)

As shown, dx1 is positive and dmu negative. Like dxmu, dx also needs a second contribution dx2, using the matrix-derivative rule derived above.

dx2 = 1. /N * np.ones((N,D)) * dmu
dx = dx1 + dx2

Combining all these steps gives the following backpropagation code.

#unfold the variables stored in cache
xhat, gamma, xmu, ivar, sqrtvar, var, eps = cache

#get the dimensions of the input/output
N,D = dout.shape

#step9
dbeta = np.sum(dout, axis=0)
dgammax = dout #not necessary, but more understandable

#step8
dgamma = np.sum(dgammax*xhat, axis=0)
dxhat = dgammax * gamma

#step7
divar = np.sum(dxhat*xmu, axis=0)
dxmu1 = dxhat * ivar

#step6
dsqrtvar = -1. /(sqrtvar**2) * divar

#step5
dvar = 0.5 * 1. /np.sqrt(var+eps) * dsqrtvar

#step4
dsq = 1. /N * np.ones((N,D)) * dvar

#step3
dxmu2 = 2 * xmu * dsq

#step2
dx1 = (dxmu1 + dxmu2)
dmu = -1 * np.sum(dxmu1+dxmu2, axis=0)

#step1
dx2 = 1. /N * np.ones((N,D)) * dmu

#step0
dx = dx1 + dx2