이전까지 했던 내용에 대한 전반적인 요약은 다음과 같다. Linear classification을 예시로 들면서 softmax, SVM(support vector machine)의 classifier를 소개했었고, 이러한 classifier가 학습되는 과정에서 사용되는 score function의 한 형태인 neural network를 언급했었다. 생체 뉴런을 유사한 형태로 표현한 perceptron의 구조와 각 연산이 가지는 의미에 대해서 살펴봤었고, 단순한 논리 구조를 벗어나 non-linearity가 적용된 여러 층의 레이어를 가지는 deep neural network를 chain-rule에 기반하여 최적화하는 과정을 통해 보다 복잡한 형태의 함수를 가지는 real-world task들에 대한 universal approximator로 사용될 수 있는 것을 확인하였다. 다만 레이어가 깊어지면 깊어질수록 representation power(복잡한 함수를 표현할 수 있는 정도)는 증가하지만 그에 따른 부작용으로 overfitting(weight값의 표준편차가 커지는 것, training dataset에 네트워크가 과적합되는 것)이 생겼고, 이를 해결하기 위한 수단으로 여러 regularization 방법들(dropout, $L_2$ regularization 등등)을 소개하였다.
이번에는 지금까지 살펴본 딥러닝 네트워크의 구조나 의의가 아닌, 실질적으로 학습 과정에서 필요한 데이터 전처리, Weight의 초기화에 대해 알아보고 학습 시 batch normalization과 같은 regularization이 가지는 장점에 대해 알아보는 글이 될 것이다.
사용할 데이터 $X$에 여러 조작을 가해보도록 하자. 여기서 사용되는 data $X$는 $X \in \mathbb{R}^{N \times D}$의 batch 단위의 matrix를 가정하도록 하자. $N$은 data sample의 갯수를 의미하고, $D$는 각 샘플의 dimensionality(차원)을 의미한다. 예를 들어 만약 $3 \times 32 \times 32$의 RGB 채널을 가지는 이미지가 총 10장이 있다면 $N = 10,~D = 3072$가 될 것이다.
가장 먼저 사용할 데이터 자체를 의미하는 original data에 대해서 살펴보자. 원본은 아무런 preprocessing이 되지 않은 그대로를 의미하고, 다음과 같이 데이터가 분포한다고 가정해보자.
간단하게 $2$차원의 데이터를 가정했으며, 좌표평면의 각 점은 sample을 의미한다.
데이터 전처리에서 사용되는 방법 중 하나는 data 분포에 존재하는 bias를 없애주기 위한 ‘mean substraction’ 작업이다. Mean substraction은 각각의 dimension에 대한 평균값을 빼주게 된다.
original_data = X #(assume that X is numpy array with shape N*D)
zero_centered_data = X - np.mean(X, axis = 0)
numpy 모듈을 잘 모르는 사람이 있을 수도 있기 때문에 언급하자면 np.mean(_, axis=0) 메소드의 경우에는 인자로 들어가게 되는 array의 $0$번째 축에 대한 평균을 구하고자 하는 것이다. $0$번째 축을 따라서 평균을 내는 것은 다시 말하자면 $2$차원의 데이터에 대해서 첫번째 column vector와 두번째 column vector($\mathbb{R}^{N \times 1}$)의 평균을 빼주는 과정이다. $2$차원 데이터에 대해서 각 column vecor의 평균으로 bias를 없애게 되면, original data가 가지던 분포의 중심점이 원점인 $(0, 0)$로 이동하게 된다.
각 차원 축은 하나의 feature로 이해할 수 있고, 만약 위의 그림과 같이 각 차원 축에 대해 분산값이 상이할 경우 학습 과정에서 최적화 시 가지는 중요도나 learning rate 비율이 달라질 수 있기 때문에 이를 어느 정도 유사하게 맞춰주는 작업이 필요하다.
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)
위에서 본 내용은 dataset sample을 정규화하는 preprocessing이었다. 일반적으로 feature vector는 위에서 보는 바와 같이 각 column마다 correlation이 어느 정도 존재하는 형태가 된다($y = ax$). 이러한 correlation을 풀어주는 작업이 PCA and Whitening이라고 생각하면 된다.
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]
연산 과정은 다음과 같다. 가지고 있는 data $X \in \mathbb{R}^{N \times D}$에 대해 zero-centered matrix를 구한 뒤, 이렇게 구해진 zero-centered matrix를 서로 inner product한 뒤 샘플의 개수 $N$으로 나눠주면 covariance matrix를 구할 수 있게 된다.
[ Cov(X) = \frac{\left(X-\mu(X) \right)^\top \left(X-\mu(X) \right)}{N}
]
구한 covariance matrix의 각 요소가 의미하는 것은 $i$번째 feature와 $j$번째 feature의 관련성을 의미한다. 위의 예시에서는 $2 \times 2$ matrix가 나오게 된다. 따라서 covariance matrix는 자동으로 symmetric matrix가 되는데, 이 행렬의 diagonal component는 feature 각각에 대한 autocorrelation이고, autocorrelation은 수학적으로 보면 variance가 된다. Covariance matrix에 singular value decomposition(SVD)를 수행하게 되면 covariance matrix로부터 eigenvector $U$, eigenvalues $V$, singular values $S$를 추출할 수 있게 된다.
U, S, V = np.linalg.svd(covariance_matrix)
SVD를 수행하는 이유는 covariance matrix의 요소가 각 feature에 대한 correlation을 표현한다고 하였는데, 이 식에서 orthonormal matrix $U$를 projection에 대한 basis로 사용하여 원본 데이터(평균에 대한 bias가 제거된)의 correlation을 없애기 위한 용도로 작용한다.
X_dr = zero_centered_data.dot(U)
Dot production이 의미하는 것은 orthonormal basis $U$에 따른 축 rotation을 의미한다. 앞서 살펴봤던 바와 같이 기존 데이터가 각 feature vector에 대한 correlation을 가지고 있었는데, 축을 회전함으로써 이를 제거해줄 수 있다(decorrelation 과정).
여기서 PCA는 Principal Component Analysis인데, 이는 SVD로 하여금 추출된 feature basis를 사용하는 것이 아니라 ‘중요한’ feature만 사용하겠다는 의미가 된다. 만약 eigenbasis에서 eigenvalue가 큰 값을 기준으로 정렬했다고 생각한다면 $D$개의 feature 중에서 정말 중요한 $100$개를 사용한다고 생각해볼 수 있다. 물론 위의 경우와 같이 $2$차원 데이터일 경우엔 굳이 PCA를 적용할 필요가 없지만, 만약 feature의 차원 수가 늘어나게 되면 유의미한 feature만 사용하는 것이 중요하다. Curse of dimension이라 표현되는 해당 문제는 data가 구성하는 manifold는 실제로 데이터가 놓인 공간이 아닌, 공간의 일부를 이루는 hyperplane 등등 특정 표면을 구성한다는 점에서 등장한다.
Xrot_reduced = np.dot(X, U[:,:100]) # Xrot_reduced becomes [N x 100]
대표적인 데이터셋 중 하나인 Swiss roll의 경우 3차원으로 구성되어있지만 데이터셋 분포는 말려있는 하나의 plane(평면)을 구성하며, 이를 그대로 3차원의 공간에 대해서 활용하는 것보다 데이터 분포를 가장 잘 나타낼 수 있는 평면을 기준으로 projection해서 사용하는 것이 PCA의 한 예시가 된다.
앞서 살펴본 예시에서 차원 수를 줄여서 사용한다면, 우리가 사용할 데이터셋은 $N \times D$에서 $N \times 100$으로 적은 feature를 가지게 되면서 그와 동시에 decorrelation이 진행된 데이터셋이 될 것이다. 여기에 추가로 data를 eigenbasis에 대해 projection 시킨 다음에 각각의 dimension을 eigenvalue로 나눠주게 되면, 각 basis에 data sample들이 가지는 variance로 normalization이 가능하다. 이를 whitening이라고 한다. $X$는 애초에 basis에 대해서 회전된 상태(projection)이고, 우리가 가지는 각 feature에 대한 분산 정보는 eigenvalue가 가지고 있다. 왜냐하면 해당 행렬의 diagonal에는 각 basis에 대한 $\sigma$가 있고, 해당 value가 covariance matrix에서 diagonal element를 결정하기 때문이다.
# whiten the data:
# divide by the eigenvalues (which are square roots of the singular values)
Xwhite = Xrot / np.sqrt(S + 1e-5)
Whitening이라고 불리는 작업은 평균이 $0$이고 identity covariance matrix를 가지는 white noise 형태의 분포를 만든다는 관점에서 나온 이름이다.
KNN(K-neareat neighbors)과 같은 instance-based learning(training dataset이 instance로 사용되어 test data의 예측에 관여하는 것)이 아닌 model-based learning은 모델 구조에 따라 weight parameter가 존재하고, 이를 update하는 과정으로 학습을 진행한다. 그렇기 때문에 학습을 시작하게 될 weight의 초기화에서도 성능의 차이가 발생하는데, 어떤 식으로 하는 것이 가장 효율적일까?
가장 먼저 생각해볼 수 있는 것은 모든 weight를 $0$으로 초기화하는 것이다. Weight의 초기화 방법은 배제하고 최적화된 네트워크의 parameter가 어떤 값을 가지는지 모두 예측할 수 없지만, 수없이 많은 parameter를 가지는 deep neural network 구조에서 weight의 절반은 positive로, 나머지는 negative로 학습이 될 것으로 예상할 수 있다(큰 수의 법칙). 그렇기에 처음부터 평균이 $0$인 상태로 시작하면 어떤 weight는 positive로, 어떤 weight는 negative로 가면서 자연스럽게 우리가 생각했던 이상적인 weight parameter의 구조를 가질 수 있지 않을까 싶지만 이는 잘못된 관점이다. 만약 network의 모든 뉴런이 같은 output을 내보내면, backpropagation 과정에서도 같은 값으로 update가 될 것이다. 즉 모든 weight가 같은 값으로 초기화되면 weight matrix가 assymetric하지 않게 된다. 간단한 예시를 위해 다음과 같은 perceptron code를 짜보았다.
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)
코드는 가독성이 떨어지기 때문에 수식으로 표현하면 다음과 같다.
[ \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}
]
위의 코드를 돌렸을 때의 weight parameter를 직접 확인해보면,
[ 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} ]
위와 같이 각 row의 모든 값이 동일하게 update되는 것을 알 수 있다. 사실상 weight parameter의 개수가 의미하는 것이 neural network에서 가용할 수 있는 node의 갯수가 되는데, 이렇게 update가 되면 굳이 큰 parameter 개수를 가지는 weight를 사용할 필요성이 없어지는 것이다. 다르게 말하면, parameter의 수가 많아도 representation power가 떨어지게 된다. 사실 이러한 문제는 모든 parameter를 $0$으로 초기화하는 상황 뿐만 아니라 같은 값으로 초기화할 때 발생하는 문제다.
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)
W1, W2 각각을 동일한 값으로 초기화한 후 학습하였다. 결과는 다음과 같다.
[ 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} ]
Overfitting 및 위에서 언급한 문제를 해결하기 위해서는 weight를 같은 값으로 초기화하는 방법을 사용할 수 없다. 다음으로 생각해볼 수 있는 상황은 모두 random하게 생성하는 것이다.
W = 0.01* np.random.randn(D, H)
여기서 $H$는 hidden layer의 output dimension을 의미한다. randn은 $D \times H$의 정규 분포 랜덤 난수를 생성하므로, 코드에 따라 $0.01$의 표준편차를 가지는 가우시안 분포를 따르는 weight matrix $W \in \mathbb{R}^{D \times H}$가 생성된다.
위의 제시된 방법은 output의 분포가 가지는 variance가 input 개수에 따라 증가할 수 있다는 것이다. 위에서와 같이 임의의 상수 0.01의 표준편차를 가지는 가우시안 분포를 생성하는 것이 아닌, 각 neuron의 output 분포를 $1$의 variance(혹은 표준편차)를 갖게끔 해준다. 특정 layer의 input 개수가 $n$이라면 해당 layer의 연산을 담당하는 neuron의 weight을
w = np.random.randn(n) / sqrt(n)
로 초기화하게 되면, 네트워크의 모든 neuron들이 같은 output distribution을 갖게되고, 이러한 방법이 convergence 속도를 증가시킬 수 있다.
[ \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}
] 해당 concept에 대한 수식 증명은 위와 같다.
위에서 언급했던 calibration 관련 문제(neuron의 output마다 분포가 달라지는 현상)를 해결하기 위해서는 모든 weight matrix를 0으로 초기화를 해야하지만, 앞서 말했던 것처럼 이대로 학습을 진행하면 모든 weight가 동일하게 학습되는 문제가 발생하기 때문에 고정된 갯수의 neuron을 정해두고 랜덤한 gaussian noise(이 때는 calibration이 진행되지 않은, 앞서 봤던 small random number에 해당된다)에서 샘플링한 weight만 연산하게 된다. 즉, 연산량을 고정시켜서 output distribution이 항상 동일하게 유지되게 하는 것이다.
Bias는 모든 값을 $0$으로 초기화해도 괜찮다. 이는 만약 weight의 assymetry가 보장되면 bias가 동일한 값을 가지더라도 node에 따른 assymetry가 유지되기 때문이다. ReLU 특성상 $0$보다 작은 값에 대해서 gradient를 주지 않기 때문에 해당 non-linearity를 가진 neuron에 대해서는 $0.01$과 같이 작은 양의 값으로 초기화하는 걸 선호하는 case도 있지만, 이러한 과정이 실질적으로 performance에 긍정적인 영향을 준다는 근거는 없으며 더 안좋은 결과를 보여주기도 한다. 따라서 bias는 $0$으로 초기화하는 것이 가장 일반적이다.
He et al. 저자들이 밝힌 바로는(참고링크) ReLU를 사용하는 network에 기반하여 다음과 같은 초기화 과정이 더 낫다고 하였다. 논문에 weight initialization 말고도 실험한 부분들이 유의미한 내용을 담고 있기 때문에 한번쯤 보는 것을 추천한다.
W = np.random.randn(n)*sqrt(2.0/n)
Deep Residual Learning for Image Recognition(이른바 ResNet)를 시작으로 convolutional neural network 구조에 batch normalization을 추가하는 것이 trend가 되었다.
Batch normalization은 말 그대로 'batch' 단위로 정규화를 진행하는 작업이며, batch normalization을 포함하여 layer normalization, instance normalization 그리고 group normalization이 있다. 각각 그림을 보게 되면 정규화를 진행하는 단위가 다른 것을 알 수 있는데, 예를 들어 layer normalization의 경우에는 단일 batch 내에서 채널 전체의 평균 및 표준편차를 통해 정규화를 진행하게 되고, Instance normalization은 단일 batch, 그리고 channel에 대해서 normalization을 진행한다. 마지막으로 group normalization은 단일 배치 내에서 채널을 묶어서 하나의 그룹을 생성하고, 이 그룹 내에서 normalization을 진행한다. 이 중에서 지금 소개할 것은 가장 흔하게 사용되는 batch normalization이다.
보통 gradient descent 알고리즘에서 단일 샘플을 사용하지 않고 mini-batch(혹은 batch라고도 부른다)를 사용하여 parameter update를 진행한다. 이유는 앞서 작성한 글에서도 확인할 수 있지만 다시 한번 언급하자면 단일 샘플로 업데이트하게 될 경우 noisy한 학습이 진행되고, batch 단위로 진행할 때와는 다르게 병렬 연산이 불가능하기 때문에 비효율적이다.
만약 여러 레이어를 가진 neural network에 batch 단위로 샘플을 통과시키게 되는 상황을 생각해보자. Batch는 랜덤하게 dataset으로부터 추출하기 때문에 각 batch 마다의 데이터 분포는 차이가 있으며, 마찬가지로 layer를 통과하면서 생기는 feature map의 분포 또한 균등하지 않다는 문제가 발생한다. 이를 ‘internal covariance shift‘라고 부른다. 사실 이 부분이 중요한 이유는 앞서 설명해왔던 weight initialization 과정과도 직결되기 때문인데, 만약 학습되는 batch 단위와 중간 layer에서의 feature map 분포 및 형태가 상이하다면 최적화 과정에서 initialization metric을 확정적으로 사용할 수 없게 된다. 따라서 이러한 batch에 따른 분포 차이를 줄여주기 위해, batch 단위로 정규화를 진행하여 변동성을 줄이고 학습을 안정화시키고자 한 것이 batch normalization의 concept이다.
각 batch에 대한 mean(평균) 및 variance(분산)은 정의에 따라 다음과 같이 구할 수 있다.
[ \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} ]
여기서 $x_i$가 의미하는 바는 $m$ 만큼의 batch가 하나의 데이터 묶음을 차지할 때, batch 내에서의 $i$번째 샘플이다. 이렇게 구한 각 batch mean, variance를 활용하여 batch dataset을 normalize 및 scale and shift를 진행한다.
[ \hat{x_i} = \frac{x_i - \mu_\text{batch}}{\sqrt{\sigma_\text{batch}^2 + \epsilon}}
]
하지만 단순히 위와 같은 방법으로 정규화를 진행한다면, non-linearity를 가지는 neural network로 하여금 linear regime에 머무르게끔 강제하는 것이 될 수 있고(representation power가 줄어드는 것을 의미), 각 레이어에서 activation function 이전의 layer의 output에 대해 적용되는 batch normalization layer이 모두 다른 형태의 representation에 대해 적용될 수 있게끔 learnable parameter인 $\gamma$, $\beta$를 두어 scale and shift를 진행한다.
[ y_i = \gamma \hat{x_i} + \beta ]
Batch normalization은 사실상 딥러닝 네트워크를 건드려본 사람이라면 필수적으로 알아야 하는 개념이지만, 말만 들어보고 적용만 했을 뿐 생각보다 실제로 구현하는 방식이나 동작하는 원리를 잘 모르고 쓰는 일이 많다. 그렇기 때문에 직접 공식을 기반으로 forward, backward가 어떻게 코드로 동작하는지 구현해보도록 하겠다. 사실 코드랑 내용 전부 다른 페이지 내용을 기반으로 작성하는 것이라 굳이 이 글을 보기 귀찮다면 참고 링크를 통해 직접 내용을 확인해봐도 좋을 것 같다.(참고링크1, 참고링크2)
배치 단위로 정규화를 진행하는 말 그대로 ‘Batch normalization’은 input $$x가 있으면 배치 단위로 정규화를 진행하는 과정이다. 이는 흔히 activation function이 뒤따르는 convolution layer 사이사이 들어가게 되며, covariance shift(feature map 및 batch 단위로 분포가 왔다 갔다 하는 현상)을 방지하여, 보다 큰 learning rate에 대해서도 안정적인 학습을 보장한다던지, 빠른 optimization을 가능하게 하는 등 다양한 장점이 있다. 하지만 Batch normalization은 무분별하게 사용되면 안되는데, 이를테면 GAN based 구조를 가진 모델에서 batch normalization이 encoder나 decoder 말단에 들어가서 학습을 방해하는 경우도 있기 때문이다. 사실 간단하게 설명하자면 batch normalization이 등장한 것은 안정적인 학습을 위함이었지만, 그렇다고 해서 모든 task에 만병통치약과 같은 존재가 될 순 없다는 것이다. 논문에 제시된 batch normalization 과정은 위에서 소개했던 수식 전개 과정과 동일하다.
[ \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} ]
가장 먼저, batch 단위로 평균($\mu$)과 분산($\sigma^2$)을 구해주고, 구한 평균과 분산을 기반으로 샘플을 정규화한다. 그런 뒤 학습 가능한 parameter인 $\gamma$와 $\beta$를 통해 scaling and shifting이 진행된다. 학습 가능한 parameter인 $\gamma$와 $\beta$는 다음과 같이 업데이트된다.
흔히 backpropagation을 계산할 때 위와 같이 diagram을 그린 뒤, chain rule에 맞춰서 local gradient를 구하고 이 값에 backward된 gradient value를 곱하게 되면, output부터 차례대로 input까지 거쳐온 모든 parameter에 대한 gradient를 구할 수 있게 된다.
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
천천히 3단계 정도씩 끊어서 볼 예정이다. 가장 먼저 샘플에 대한 평균을 구한다. sample x의 모양은 batch * dimension이므로 axis=0을 기준으로 평균을 구해야 batch에 대한 평균을 구할 수 있다. 참고로 np.sum() 메소드가 아닌 np.mean() 메소드를 통해 평균값을 바로 구할 수 있다.
mu = np.mean(x, axis=0)
평균을 구했으니 이를 sample에서 빼준다. 여기서 중요한 점은 우리가 x - mu를 두 가지 계산에 동시에 활용할 것인데, 바로 첫 번째는 아래와 같이 normalize된 $x$를 구할 때의 x - mu term에 활용할 것이고 동시에 var를 구하기 위한 식에도 적용할 것이다.
[ \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} ]
위의 두 식에 공통적으로 들어있는 $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는 위에서 구했던 x - mu의 제곱이기 때문에 이를 평균낸 것이 분산 공식이다. 위에서 언급했던 것과 마찬가지로 이 식도 np.mean() 메소드로 대체 가능하다.
var = np.mean(sq, axis = 0)
그리고 step 5에서 eps를 더하는 형태의 방법은 흔히 컴퓨팅 환경에서 0으로 division되는 error를 방지하기 위함이다. 나머지 계산을 통해 구하고자 하는 output을 나타내면 다음과 같다.
#step7: execute normalization
xhat = xmu * ivar
#step8: Nor the two transformation steps
gammax = gamma * xhat
#step9
out = gammax + beta
xhat(normalized된 x), 학습 가능한 parameter인 gamma, beta를 통한 scaling과 shifting 과정이다.
[ \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} ]
사실 backward propation 과정은 위와 같이 논문에 수식으로 정리되어 있다. 하지만 실제로 이를 코드로 구현하는 작업이 생각보다는 복잡하다. 각각의 식이 유도된 부분을 각 gate마다 천천히 살펴보면 다음과 같다. 우선 가장 말단에 있는 gate부터 backpropagation을 진행하게 되면,
[ \text{out} = \gamma \hat{x} + \beta
] 위와 같은 연산을 하게 되므로 덧셈 연산 gate의 두 input인 $\gamma \hat{x},~\beta$에 대한 local gradient는 각각 다음과 같이 구할 수 있다.
dgammax = dout
dbeta = np.sum(dout, axis=0)
dbeta를 구하는 과정에서 np.sum() 메소드가 사용되는 부분이 이해가 잘 가지 않아서 살펴보니, numpy broadcasting은 사실 $\beta$의 크기 그대로를 더하는 것이 아니라, 이 값을 broadcasting하여 gammax의 dimension에 맞게끔 확장시켜 더하게 된다. 따라서 dbeta를 구하는 것은 broadcasting을 고려해야하므로 이처럼 표현된 것으로 이해할 수 있었다.
그 다음은 곱셈 연산 게이트에 대해서 backpropagation을 진행한다. 이 operation에 대한 backpropagation을 계산하기 전, dgammax가 이전의 backpropagation value이기 때문에 chain rule에 따라, 해당 value에 local gradient를 곱하여 해당 gate의 각 input에 대해 흘러가는 gradient를 구하게 된다. 사실 * operation에 대한 gate local gradient는 이전 게시글에서 다룬 것과 같이, gate로 들어오는 backpropagation value에 각자 다른 길로 들어오는 input value를 곱해주면 된다. 이를 테면 위쪽 input은 xhat이고 아래쪽 input은 gamma이므로, 위쪽 gate에 대한 backpropagation value는 dgammax에 gamma를 곱한 결과가 dxhat이 되고 반대로 아래쪽 gate에 대한 backpropagation value는 dgammax에 xhat을 곱한 결과가 된다.
dxhat = dgammax * gamma
dgamma = np.sum(dgammax * xhat, axis=0)
그리고 앞서 살펴본 바와 같이 gamma 또한 beta와 같이 broadcasting되는 값이기 때문에 np.sum() 메소드를 통해 dimension을 맞춰주게 된다. 이제 학습 가능한 parameter인 gamma와 beta에 대한 backpropagation이 끝났고, normalization에서의 backpropagation을 쭉 계산해보면 다음과 같다.
가장 오른쪽의 곱셈 연산 gate에 대해서는 앞서 했던 연산과 동일하므로 주어진 코드에서 연산은 다음과 같이 진행된다.
divar = np.sum(dxhat*xmu, axis=0)
dxmu1 = dxhat * ivar
sample variance divar(inverse variance) 또한 broadcasting 되었으므로 np.sum 해주고, dxmu1은 이전의 backpropagation value였던 dxhat에 ivar를 곱해준 것과 같다. 여기서 중요한 것은 dxmu에 대한 계산이 두 개가 필요하다는 것인데, 이는 아래와 같이 x - mu가 두 갈래로 나뉘어 계산되었기 때문이다.
따라서 이 gate에 대한 backprop을 계산하기 위해서는 분리된 두 input에 대한 gradient를 따로 계산한 뒤에 더해줘야 한다. 위로 가는 dxmu1은 계산했고 나머지 아래 부분에 대한 공식을 차례로 보면, inverse에 대한 local gradient는 $\frac{d}{dx} \left( \frac{1}{x} \right) = -\frac{1}{x^2}$이기 때문에,
dsqrtvar = -1. /(sqrtvar**2) * divar
코드가 위와 같으며 마찬가지로 square root에 대한 local gradient는 $\frac{d}{dx} \left( x + \epsilon \right) = \frac{1}{2 \sqrt{x + \epsilon}}$이 되기 때문에,
dvar = 0.5 * 1. /np.sqrt(var+eps) * dsqrtvar
이처럼 표현 가능하다. $\frac{1}{N} \sum_i x_i$에 대한 부분이 조금 복잡한데, 현재 미분의 대상이 되는 것이 행렬이기 때문에 다음과 같이 정의할 수 있다.
[ \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
마지막 부분은 $\frac{d}{dx} \left( x^2 \right) = 2x$에서,
dxmu2 = 2 * xmu * dsq
이제 드디어 dxmu2를 구했으므로 뺄셈 연산 gate에 대한 local gradient를 구하게 되면 mu는 빼주고 x는 더해주는 process이므로,
dx1 = (dxmu1 + dxmu2)
dmu = -1 * np.sum(dxmu1 + dxmu2, axis=0)
위에서 보는 바와 같이 dx1은 positive(+), dmu는 negative(-) 방향이 된다. dx 또한 dxmu 계산과 동일하게 dx2 연산이 추가로 필요한데, 이는 앞서 구했던 matrix의 미분 공식과 같은 공식이 적용된다.
dx2 = 1. /N * np.ones((N,D)) * dmu
dx = dx1 + dx2
이를 모두 합친 과정이 다음과 같은 backpropagation 코드가 된다.
#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