Forma matricial de retropropagación con normalización de lotes

La normalización por lotes se ha acreditado con mejoras sustanciales de rendimiento en redes neuronales profundas. Un montón de material en Internet muestra cómo implementarlo en una activación por activación. Ya implementé backprop usando álgebra matricial, y dado que estoy trabajando en lenguajes de alto nivel (mientras confío en Rcpp(y eventualmente GPU) para la multiplicación de matriz densa), arrancar todo y recurrir a for-loops probablemente ralentizaría mi código sustancialmente, además de ser un gran dolor.

La función de normalización por lotes es

$$b(x_p) = \gamma \left(x_p - \mu_{x_p}\right) \sigma^{-1}_{x_p} + \beta$$ donde

• $x_p$ es el$p$ th nodo, antes de que se active
• $\gamma$ y$\beta$ son parámetros escalares
• $\mu_{x_p}$ y$\sigma_{x_p}$ son la media y la DE de$x_p$ . (Tenga en cuenta que la raíz cuadrada de la varianza más un factor de fudge se usa normalmente; supongamos que los elementos no son cero para la compacidad)

En forma de matriz, la normalización de lotes para una capa completa sería donde

$$b(\mathbf{X}) = \left(\gamma\otimes\mathbf{1}_p\right)\odot \left(\mathbf{X} - \mu_{\mathbf{X}}\right) \odot\sigma^{-1}_{\mathbf{X}} + \left(\beta\otimes\mathbf{1}_p\right)$$

• es N × $\mathbf{X}$$N\times p$
• es un vector columna de unos$\mathbf{1}_N$
• y β son ahora fila p- vectores de los parámetros de normalización por capa$\gamma$$\beta$$p$
• y σ X sonmatrices N × p , donde cada columna es unvector N de medias en columna y desviaciones estándar$\mu_{\mathbf{X}}$$\sigma_{\mathbf{X}}$$N \times p$$N$
• es el producto Kronecker y ⊙ es el producto elementwise (Hadamard)$\otimes$$\odot$

Una red neuronal de una capa muy simple sin normalización por lotes y un resultado continuo es

$$y = a\left(\mathbf{X\Gamma}_1\right)\Gamma_2 + \epsilon$$

dónde

• es p 1 × p $\Gamma_1$$p_1 \times p_2$
• es p 2 × $\Gamma_2$$p_2 \times 1$
• es la función de activación$a(.)$

Si la pérdida es , a continuación, los gradientes son ∂ $R = N^{-1}\displaystyle\sum\left(y - \hat{y}\right)^2$

$$\begin{array}{lr} \frac{\partial R}{\partial \Gamma_1} = -2\mathbf{V}^T \hat\epsilon\\ \frac{\partial R}{\partial \Gamma_2} = \mathbf{X}^T \left(a'(\mathbf{X}\mathbf{\Gamma}_1) \odot -2\hat\epsilon \mathbf{\Gamma}_2^T\right) \\ \end{array}$$

dónde

• $\mathbf{V} = a\left(\mathbf{X}\Gamma_1\right)$
• $\hat{\epsilon} = y-\hat{y}$

Bajo la normalización de lotes, la red se convierte en o

$$y = a\left(b\left(\mathbf{X}\Gamma_1\right)\right)\Gamma_2$$ $$y = a\Big(\left(\gamma\otimes\mathbf{1}_N\right)\odot \left(\mathbf{X\Gamma_1} - \mu_{\mathbf{X\Gamma_1}}\right) \odot\sigma^{-1}_{\mathbf{X\Gamma_1}} + \left(\beta\otimes\mathbf{1}_N\right)\Big)\mathbf{\Gamma_2}$$ No tengo idea de cómo calcular los derivados de los productos Hadamard y Kronecker. Sobre el tema de los productos Kronecker, la literatura se vuelve bastante misteriosa.

¿Existe una forma práctica de calcular , ∂ R / ∂ β y ∂ R / ∂ Γ 1 dentro del marco de la matriz? ¿Una expresión simple, sin recurrir a la computación nodo por nodo?$\partial R/\partial \gamma$$\partial R/\partial \beta$$\partial R/\partial \mathbf{\Gamma_1}$

Actualización 1:

He descubierto - más o menos. Es: 1 T N ( un ' ( X Γ 1 ) ⊙ - 2 varepsilon Γ T 2 ) Algunos R código demuestra que esto es equivalente a la forma de bucle para hacerlo. Primero configure los datos falsos:$\partial R/\partial \beta$

$$\mathbf{1}_{N}^T \left(a'(\mathbf{X}\mathbf{\Gamma}_1) \odot -2\hat\epsilon \mathbf{\Gamma}_2^T\right)$$

set.seed(1)
library(dplyr)
library(foreach)

#numbers of obs, variables, and hidden layers
N <- 10
p1 <- 7
p2 <- 4
a <- function (v) {
  v[v < 0] <- 0
  v
}
ap <- function (v) {
  v[v < 0] <- 0
  v[v >= 0] <- 1
  v
}

# parameters
G1 <- matrix(rnorm(p1*p2), nrow = p1)
G2 <- rnorm(p2)
gamma <- 1:p2+1
beta <- (1:p2+1)*-1
# error
u <- rnorm(10)

# matrix batch norm function
b <- function(x, bet = beta, gam = gamma){
  xs <- scale(x)
  gk <- t(matrix(gam)) %x% matrix(rep(1, N))
  bk <- t(matrix(bet)) %x% matrix(rep(1, N))
  gk*xs+bk
}
# activation-wise batch norm function
bi <- function(x, i){
  xs <- scale(x)
  gk <- t(matrix(gamma[i]))
  bk <- t(matrix(beta[i]))
  suppressWarnings(gk*xs[,i]+bk)
}

X <- round(runif(N*p1, -5, 5)) %>% matrix(nrow = N)
# the neural net
y <- a(b(X %*% G1)) %*% G2 + u

Luego calcule derivados:

# drdbeta -- the matrix way
drdb <- matrix(rep(1, N*1), nrow = 1) %*% (-2*u %*% t(G2) * ap(b(X%*%G1)))
drdb
           [,1]      [,2]    [,3]        [,4]
[1,] -0.4460901 0.3899186 1.26758 -0.09589582
# the looping way
foreach(i = 1:4, .combine = c) %do%{
  sum(-2*u*matrix(ap(bi(X[,i, drop = FALSE]%*%G1[i,], i)))*G2[i])
}
[1] -0.44609015  0.38991862  1.26758024 -0.09589582

Ellos coinciden. Pero todavía estoy confundido, porque realmente no sé por qué esto funciona. Las notas de MatCalc a las que hace referencia @Mark L. Stone dicen que la derivada de debe ser$\beta \otimes \mathbf{1}_N$

$$\frac{\partial A \otimes B}{\partial A} = \left(I_{nq} \otimes T_{mp}\right)\left(I_n\otimes vec(B) \otimes I_m\right)$$ $m$$n$$p$$q$$A$$B$$T$

# playing with the kroneker derivative rule
A <- t(matrix(beta)) 
B <- matrix(rep(1, N))
diag(rep(1, ncol(A) *ncol(B))) %*% diag(rep(1, ncol(A))) %x% (B) %x% diag(nrow(A))
     [,1] [,2] [,3] [,4]
 [1,]    1    0    0    0
 [2,]    1    0    0    0
 snip
[13,]    0    1    0    0
[14,]    0    1    0    0
snip
[28,]    0    0    1    0
[29,]    0    0    1    0
[snip
[39,]    0    0    0    1
[40,]    0    0    0    1

$\gamma$$\mathbf{\Gamma_1}$$\beta \otimes \mathbf{1}$

Actualización 2

Leyendo libros de texto, estoy bastante seguro de que $\partial R/\partial \Gamma_1$$\partial R/\partial \gamma$vec()$\partial R/\partial \Gamma_1$$w\odot\mathbf{X\Gamma_1}$$\mathbf{\Gamma_1}$$w \equiv (\gamma \otimes \mathbf{1}) \odot \sigma_{\mathbf{X\Gamma_1}}^{-1}$

$w\odot\mathbf{X}$$w$$\mathbf{X}$

$$\partial(A \odot B) = \partial A \odot B + A \odot \partial B$$

y de esto , eso

$$\frac{\partial vec(w \odot \mathbf{X\Gamma_1})}{\partial vec(\mathbf{\Gamma_1})^T} = vec(\mathbf{X\Gamma_1})I\frac{\partial vec(w)}{\partial vec(\mathbf{\Gamma_1})^T} + vec(w)I\frac{\partial vec(\mathbf{X\Gamma_1})}{\partial vec(\mathbf{\Gamma_1})^T}$$

Actualización 3

Progresando aquí. Me desperté a las 2 de la madrugada con esta idea. Las matemáticas no son buenas para dormir.

$\partial R/\partial \mathbf{\Gamma_1}$

• $w \equiv (\gamma \otimes \mathbf{1}) \odot \sigma_{\mathbf{X\Gamma_1}}^{-1}$
• $\text{"stub"} \equiv a'(b(\mathbf{X\Gamma}_1)) \odot -2\hat\epsilon \mathbf{\Gamma}_2^T$

$$\frac{\partial R}{\partial \Gamma_1} = \frac{\partial w \odot \mathbf{X\Gamma}_1}{\partial \Gamma_1}\left(\text{"stub"}\right)$$ $i$$j$$\mathbf{I}$ $$\frac{\partial R}{\partial \Gamma_{ij}} = \left(w_i \odot \mathbf{X_i}\right)^T\left(\text{"stub"}_j\right)$$ $$\frac{\partial R}{\partial \Gamma_{ij}} = \left(\mathbf{I} w_i \mathbf{X_i}\right)^T\left(\text{"stub"}_j\right)$$ $$\frac{\partial R}{\partial \Gamma_{ij}} = \mathbf{X_i}^T\mathbf{I} w_i\left(\text{"stub"}_j\right)$$ $$\frac{\partial R}{\partial \Gamma} = \mathbf{X}^T\left(\text{"stub"}\odot w\right)$$

Y, de hecho, es:

stub <- (-2*u %*% t(G2) * ap(b(X%*%G1)))
w <- t(matrix(gamma)) %x% matrix(rep(1, N)) * (apply(X%*%G1, 2, sd) %>% t %x% matrix(rep(1, N)))
drdG1 <- t(X) %*% (stub*w)

loop_drdG1 <- drdG1*NA
for (i in 1:7){
  for (j in 1:4){
    loop_drdG1[i,j] <- t(X[,i]) %*% diag(w[,j]) %*% (stub[,j])
  }
}

> loop_drdG1
           [,1]       [,2]       [,3]       [,4]
[1,] -61.531877  122.66157  360.08132 -51.666215
[2,]   7.047767  -14.04947  -41.24316   5.917769
[3,] 124.157678 -247.50384 -726.56422 104.250961
[4,]  44.151682  -88.01478 -258.37333  37.072659
[5,]  22.478082  -44.80924 -131.54056  18.874078
[6,]  22.098857  -44.05327 -129.32135  18.555655
[7,]  79.617345 -158.71430 -465.91653  66.851965
> drdG1
           [,1]       [,2]       [,3]       [,4]
[1,] -61.531877  122.66157  360.08132 -51.666215
[2,]   7.047767  -14.04947  -41.24316   5.917769
[3,] 124.157678 -247.50384 -726.56422 104.250961
[4,]  44.151682  -88.01478 -258.37333  37.072659
[5,]  22.478082  -44.80924 -131.54056  18.874078
[6,]  22.098857  -44.05327 -129.32135  18.555655
[7,]  79.617345 -158.71430 -465.91653  66.851965

Actualización 4

$\partial R / \partial \gamma$

• $\widetilde{\mathbf{X\Gamma}} \equiv \left(\mathbf{X\Gamma} - \mu_{\mathbf{X\Gamma}}\right)\odot \sigma^{-1}_\mathbf{X\Gamma}$
• $\tilde\gamma \equiv \gamma \otimes\mathbf{1}_N$

$$\frac{\partial R}{\partial \tilde\gamma} = \frac{\partial \tilde\gamma \odot \widetilde{\mathbf{X\Gamma}}}{\partial \tilde\gamma}\left(\text{"stub"}\right)$$ $$\frac{\partial R}{\partial \tilde\gamma_i} = (\widetilde{\mathbf{X\Gamma}})_i^T \mathbf{I}\tilde\gamma_i \left(\text{"stub"}_i\right)$$ $$\frac{\partial R}{\partial \tilde\gamma} = (\widetilde{\mathbf{X\Gamma}})^T \left(\text{"stub"} \odot \tilde\gamma \right)$$

Más o menos coincide:

drdg <- t(scale(X %*% G1)) %*% (stub * t(matrix(gamma)) %x% matrix(rep(1, N)))

loop_drdg <- foreach(i = 1:4, .combine = c) %do% {
  t(scale(X %*% G1)[,i]) %*% (stub[,i, drop = F] * gamma[i])  
}

> drdg
           [,1]      [,2]       [,3]       [,4]
[1,]  0.8580574 -1.125017  -4.876398  0.4611406
[2,] -4.5463304  5.960787  25.837103 -2.4433071
[3,]  2.0706860 -2.714919 -11.767849  1.1128364
[4,] -8.5641868 11.228681  48.670853 -4.6025996
> loop_drdg
[1]   0.8580574   5.9607870 -11.7678486  -4.6025996

$\gamma$

Parece que he respondido a mi propia pregunta, pero no estoy seguro de si estoy en lo correcto. En este punto, aceptaré una respuesta que prueba rigurosamente (o refuta) lo que he pirateado juntos.

while(not_answered){
  print("Bueller?")
  Sys.sleep(1)
}

No es una respuesta completa, pero para demostrar lo que sugerí en mi comentario si 2 , ...

$$b(X)=(X−e_N\mu_X^T)ΓΣ_X^{-1/2}+e_N\beta^T$$ $\Gamma=\mathop{\mathrm{diag}}(\gamma)$$\Sigma_X^{-1/2}=\mathop{\mathrm{diag}}(\sigma_{X_1}^{-1},\sigma_{X_2}^{-1},\dots)$$e_N$ $$\nabla_\beta R=[-2\hat{\epsilon}(\Gamma_2^T\otimes I)J_X(a)(I\otimes e_N)]^T$$ $-2\hat{\epsilon}(\Gamma_2^T\otimes I)=\mathop{\mathrm{vec}}(-2\hat{\epsilon}\Gamma_2^T)^T$$J_X(a)=\mathop{\mathrm{diag}}(\mathop{\mathrm{vec}}(a^\prime(b(X\Gamma_1))))$ $$\nabla_\beta R=(I\otimes e_N^T)\mathop{\mathrm{vec}}(a^\prime(b(X\Gamma_1))\odot-2\hat{\epsilon}\Gamma_2^T)=e_N^T(a^\prime(b(X\Gamma_1))\odot-2\hat{\epsilon}\Gamma_2^T)$$ $\mathop{\mathrm{vec}}(AXB)=(B^T\otimes A)\mathop{\mathrm{vec}}(X)$ $$\begin{align}\nabla_\gamma R&=[-2\hat{\epsilon}(\Gamma_2^T\otimes I)J_X(a)(\Sigma_{X\Gamma_1}^{-1/2}\otimes (X\Gamma_1-e_N\mu_{X\Gamma_1}^T))K]^T\\&=K^T\mathop{\mathrm{vec}}((X\Gamma_1-e_N\mu_{X\Gamma_1}^T)^TW\Sigma^{-1/2}_{X\Gamma_1})\\&=\mathop{\mathrm{diag}}((X\Gamma_1-e_N\mu_{X\Gamma_1}^T)^TW\Sigma^{-1/2}_{X\Gamma_1})\end{align}$$ $W=a^\prime(b(X\Gamma_1))\odot-2\hat{\epsilon}\Gamma_2^T$$K$$Np\times p$$d\Gamma_{i\neq j}=0$$b$$\gamma_i$$\gamma_i$$\Gamma_1$$w$$\Sigma_X$$\mu_X$$X$