Restricted Boltzmann Machines (RBM) con Tensorflow

Hoy vamos a hablar de unas redes muy distintas a las que habitualmente hablamos, estas son las Restricted Boltzmann Machines (RBM).
Las RBM son redes neuronales estocásticas generativas y habitualmente se usan tanto para problemas supervisados, semi-supervisados o no supervisados. Cuando decimos generativa se refiere a que estima la distribución poblacional mediante las muestras y cuando decimos estocástica indica que dentro de la red se introducen variables aleatorias.
¿Porque he subrayado semi-supervisados  o no supervisados? Porque resolver problemas supervisados es bastante «sencillo», pero cuando no tenemos datos etiquetados o los tenemos parcialmente etiquetados las cosas se nos puede complicar bastante.
Este tipo de redes nos permiten tratar este tipo de aprendizaje. Serán capaces de aprender características sobre datos que no tenemos etiqueta alguna o en el que solo tenemos un pequeño conjunto de muestras etiquetados.

Estructura de la red

Se conoce como visible units las unidades de la red que sí conocemos y hidden units las unidades que desconocemos o tenemos parcialmente etiquetadas.

Si os fijáis en la anterior figura no hay relaciones entre las visible units (entre ellas) o entre las hidden units (entre ellas), esta restricción es la que diferencia las Restricted Boltzmann Machines de las Boltzmann Machines convencionales (intratables por el coste computacional que suponen).
Las Restricted Boltzmann Machines nos permite asumir que cada visible unit es independiente las demás y que lo mismo sucede para las hidden units.
Este tipo de red proviene conceptualmente de la física, donde existe una función que calcula la energía total de un sistema.
E(v, h) = - \sum_i a_i v_i - \sum_j b_j h_j - \sum_{i,j} v_i W_{i,j} h_j Donde v serán los datos que sí conocemos (visible), h (hidden) serán la cantidad de características que queremos obtener o las características parcialmente etiquetadas.
En los parámetros que encontraremos són:

• a_i regula el bias, importancia o contribución de la v_i en la energía total del sistema.
• b_j regula el bias, importancia o contribución de la h_j en la energía total del sistema.
• W_{i,j} regula la relación entre las v y h.

Podemos estimar la probabilidad conjunta de las visible units/hidden units a partir de la energía anteriormente calculada:
P(v, h) = \frac{1}{Z} e^{- E(v, h)} Z, es un factor que normaliza para que el valor sea de [0, 1].
Z= \sum_{v, h} e^{- E(v, h)} Cómo podéis ver [/latex]Z es extremadamente complicada de estimar, supone probar todos los posibles valores tanto de [latex]v como de h, si tenemos 100 visible units y 100 hidden units binarias hablamos de 2^{100 + 100} combinaciones distintas. Lo que produce una complejidad exponencial y es muy difícil de tratar directamente...
La probabilidad de las visible units tiene exactamente el mismo problema:
P(v) = \frac{1}{Z} \sum_h e^{- E(v, h)}

Estimar la probabilidad condicionada

¿Como estimamos las probabilidades condicionadas? Centrémonos en P(v | h):
P(h | v) =\frac{\cancel{\frac{1}{Z}} e^{- E(v, h)}}{\cancel{\frac{1}{Z}} \sum_h e^{- E(v, h)}} podemos cancelar a \frac{1}{Z}. Y desarrollar E(v, h):
= \frac{e^{- E(v, h)}}{\sum_h e^{- E(v, h)}} = \frac{e^{\sum_i a_i v_i + \sum_j b_j h_j + \sum_{i,j} v_i W_{i,j} h_j}}{\sum_h e^{\sum_i a_i v_i + \sum_j b_j h_j + \sum_{i,j} v_i W_{i,j} h_j}} podemos cancelar el termino e^{\sum_i a_i v_i} ya que no depende de h.
= \frac{\cancel{e^{\sum_i a_i v_i}} e^{\sum_j b_j h_j + \sum_{i,j} v_i W_{i,j} h_j}}{\cancel{e^{\sum_i a_i v_i}}  \sum_h e^{\sum_j b_j h_j + \sum_{i,j} v_i W_{i,j} h_j}} podemos convertir los sumatorios dentro de la e por productorios:
=\frac{ \prod_j e^{b_j h_j + \sum_{i} v_i W_{i,j} h_j}}{\sum_h \prod_j e^{b_j h_j + \sum_{i} v_i W_{i,j} h_j}} como el sumatorio del producto esta en términos de los indices y son independientes entre sí, lo podemos convertir al producto de un sumatorio:
= \frac{ \prod_j e^{b_j h_j + \sum_{i} v_i W_{i,j} h_j}}{\sum_{h_j=0} \dots \sum_{h_j=n} (e^{b_0 h_0 + \sum_{i} v_i W_{i,0} h_0}) \dots (e^{b_n h_n + \sum_{i} v_i W_{i,n} h_n})} = \frac{ \prod_j e^{b_j h_j + \sum_{i} v_i W_{i,j} h_j}}{\sum_{h_j=0} (e^{b_0 h_0 + \sum_{i} v_i W_{i,0} h_0}) \dots \sum_{h_j=n} (e^{b_n h_n + \sum_{i} v_i W_{i,n} h_n})} fijaros como cada elemento del productorio solo depende de un solo valor h_j del sumatorio, esto nos permite re-ordenar:
= \frac{ \prod_j e^{b_j h_j + \sum_{i} v_i W_{i,j} h_j}}{ \prod_j \sum_h e^{b_j h_j + \sum_{i} v_i W_{i,j} h_j}} esto nos permite decir que P(h_j | v) es:
P(h | v) = \prod_j \underbrace{\frac{e^{b_j h_j + \sum_{i} v_i W_{i,j} h_j}}{\sum_{h_j} e^{b_j h_j + \sum_{i} v_i W_{i,j} h_j}}}_{P(h_j | v)} P(v | h) =\frac{\cancel{\frac{1}{Z}} e^{- E(v, h)}}{\cancel{\frac{1}{Z}} \sum_{v} e^{- E(v, h)}}=\frac{e^{- E(v, h)}}{\sum_v e^{- E(v, h)}}= \prod_i \underbrace{\frac{e^{a_i v_i + \sum_{j} v_i W_{i,j} h_j}}{\sum_{v_i} e^{a_i v_i + \sum_{j} v_i W_{i,j} h_j}}}_{P(v_i | h)}

Otra forma de enfocarlo es partir de E(v, h), desarrollarlo para el caso de una sola hidden unit, calcular la probabilidad condicionada y obtendremos:
P(h_{unica} | v) = \frac{e^{b_{unica} h_{unica} + \sum_{i}  v_i W_i h_{unica}}}{\sum_h e^{b_{unica} h_{unica} + \sum_{i}  v_i W_i h_{unica}}} Después teniendo en consideración que en la formula E(v, h), las unidades h entre sí son independientes, \forall_{i,j} | i \neq j, P(h_i, h_j | v) = P(h_i | v) \cdot P(h_j | v). Llegaríamos a la misma conclusión.

Considerando las unidades valores binarios

Cuando tenemos en cuenta que las hidden units solo pueden comprender valores discretos binarios [0, 1], las cosas se nos simplifican considerablemente (centrémonos en P(h_j = 1 | v)):
P(h_j = 1 | v) =\frac{e^{b_j h_j + \sum_{i} v_i W_{i,j} h_j}}{\sum_{h_j} e^{b_j h_j + \sum_{i} v_i W_{i,j} h_j}} =\frac{e^{b_j + \sum_{i} v_i W_{i,j}}}{\underbrace{ e^{b_j 0 + \sum_{i} v_i W_{i,j}0}}_ { 1} + e^{b_j 1 + \sum_{i} v_i W_{i,j} 1}} =\frac{e^{b_j + \sum_{i} v_i W_{i,j}}}{1 +e^{b_j + \sum_{i} v_i W_{i,j}}} P(h_j = 1 | v) = sigmoid \left ( b_j + \sum_{i} v_i W_{i,j} \right ) P(v_i = 1 | h) = sigmoid \left ( a_i + \sum_{j} W_{i,j} h_j \right )

Considerando las unidades visibles K valores discretos

Esta vez los valores de las visible units pueden comprender K distintos valores discretos, tranformando RBM en un problema multinomial. Los K valores estarán codificados usando one-hot-encode.
Para ello ahora deberemos considerar que para un valor cualquiera de i y j los parámetros a_j y W_{i,j} ya no serán un simple valor sino que serán un vector de K posiciones. Entonces:

• a: es una matriz num(visibles) \times K.
• W: es un cubo num(visibles) \times num(hiddens) \times K.

P(h_j = k | h) =\frac{e^{b_{j, k} + \sum_{i} v_i W_{i,j,k}}}{\sum_k e^{b_{j, k} + \sum_{i} v_i W_{i,j,k}}} P(v_i = k | h) = P(v_{i,k} = 1 | h) = softmax_i \left( {a_{i, k} + \sum_{j} h_i W_{i,j,k}} \right ) P(h_j = 1 | v) = sigmoid \left( {b_j + \sum_{i} \sum_k  v_{i, k} W_{i,j,k}} \right )

Aprendiendo los pesos

Si habéis llegado hasta aquí pensareis: ¿Muy bien, la formula es muy bonita pero como puede aprender por si sola los pesos de h? Resolver RBM implicará buscar los parámetros que maximicen P(v | a, b, W), estaremos diciendo a nuestra función de energía: "Eh por favor, buscame los parámetros en los que v y h pierdan la menos información posible" y como consecuencia encontraremos una representación más compacta en h de v. Como siempre, no usaremos una única muestra para realizar la estimación, sino que usaremos N, quedando la solución de nuestro problema de la siguiente forma (es un log-MLE):
argmax_{a, b, W} \sum_i^N log P(v^i | a, b, W) Si desglosamos un poco la formula obtenemos la probabilidad conjunta:
argmax_{a, b, W} \sum_i^N log \sum_h P(v^i, h | a, b, W) argmax_{a, b, W} \sum_i^N log \sum_h\frac{1}{Z} e^{- E(v^i, h)} como Z no depende de h y es constante para todas las N muestras...
argmax_{a, b, W} \sum_i^N log \sum_h e^{- E(v^i, h)} - N log Z argmax_{a, b, W} \underbrace{ \sum_i^N log \sum_h e^{- E(v^i, h)}}_{\text{fase positiva}} - \underbrace{ N log \sum_{v, h} e^{- E(v, h)}}_{\text{fase negativa}} La fase positiva aumenta la probabilidad de los h creados por los ejemplos y la fase negativa disminuye la probabilidad de aquellos h que han sido creados "artificialmente". Con esto se trata de eliminar el ruido que genera propiamente la RBM.

Derivando

Ahora toca el paso más engorroso de todos (más aún) vamos a derivar respecto los parámetros para encontrar como corregirlos mediante el descenso de gradiente.

Fase positiva

Respecto a
\frac{\partial \sum_i^N log \sum_h e^{- E(v^i, h)}}{\partial a} = \sum_i^N \underbrace{\frac{1}{\sum_h e^{- E(v^i, h)}} \sum_h \left [ e^{- E(v^i, h)} \underbrace{\frac{- \partial E(v^i, h)}{\partial a}}_{v^i} \right ]}_{v^i} Si analizamos un poco la formula vemos que se está haciendo la esperanza de v^i respecto la probabilidad de todos los valores distintos que puede tener h sujeto a v, es decir: E_{P(h| v^i)}[v^i]. Y v^i es constante para todos estos valores. De hecho v no depende de \sum h de modo que podemos sacar v fuera del sumatorio.
\frac{\partial \sum_i^N log \sum_h e^{- E(v^i, h)}}{\partial a} = \sum_i^N {v^i \underbrace{\cancel{\frac{1}{\sum_h e^{- E(v^i, h)}} \sum_h e^{- E(v^i, h)}}}_{1}} \frac{\partial \sum_i^N log \sum_h e^{- E(v^i, h)}}{\partial a} = \sum_i^N v^i Respecto b
\frac{\partial \sum_i^N log \sum_h e^{- E(v^i, h)}}{\partial b} = \sum_i^N \underbrace{\frac{1}{\sum_h e^{- E(v^i, h)}} \sum_h \left [ e^{- E(v^i, h)} \underbrace{\frac{- \partial E(v^i, h)}{\partial b}}_{h} \right ]}_{E_{P(h| v^i)}[h]} Podemos pensar... ¡Que narices es E_{P(h| v^i)}[h]! Recordar que los valores de h son siempre binarios de modo que E_{P(h| v^i)}[h] = P(h| v^i).
\frac{\partial \sum_i^N log \sum_h e^{- E(v^i, h)}}{\partial b} = P(h| v^i) Respecto W
\frac{\partial \sum_i^N log \sum_h e^{- E(v^i, h)}}{\partial W} = \sum_i^N \underbrace{\frac{1}{\sum_h e^{- E(v^i, h)}} \sum_h \left [ e^{- E(v^i, h)} \underbrace{\frac{- \partial E(v^i, h)}{\partial W}}_{v^i h} \right ]}_{v^i P(h| v^i)} \frac{\partial \sum_i^N log \sum_h e^{- E(v^i, h)}}{\partial W} = v^i P(h| v^i)

Fase negativa

Respecto a
\frac{\partial N log \sum_{v, h} e^{- E(v, h)}}{\partial a} = N \underbrace{\frac{1}{\sum_{v, h} e^{- E(v, h)}}\sum_{v, h} \left [e^{- E(v, h)} \underbrace{\frac{- \partial E(v, h)}{\partial a}}_{v} \right ]}_{E_{P(v, h)}[v]} \frac{\partial N log \sum_{v, h} e^{- E(v, h)}}{\partial a} = N E_{P(v, h)}[v] Respecto b
\frac{\partial N log \sum_{v, h} e^{- E(v, h)}}{\partial b} = N \underbrace{\frac{1}{\sum_{v, h} e^{- E(v, h)}}\sum_{v, h} \left [e^{- E(v, h)} \underbrace{\frac{- \partial E(v, h)}{\partial b}}_{v} \right ]}_{E_{P(v, h)}[h]} \frac{\partial N log \sum_{v, h} e^{- E(v, h)}}{\partial b} = N E_{P(v^i, h)}[h] Respecto W
\frac{\partial N log \sum_{v, h} e^{- E(v, h)}}{\partial W} = N \underbrace{\frac{1}{\sum_{v, h} e^{- E(v, h)}}\sum_{v, h} \left [e^{- E(v, h)} \underbrace{\frac{- \partial E(v, h)}{\partial W}}_{v h} \right ]}_{E_{P(v, h)}[v^T h]} \frac{\partial N log \sum_{v, h} e^{- E(v, h)}}{\partial W} = N E_{P(v^i, h)}[v^T h]

Juntamos ambas partes

\nabla_a = \sum_i^N v^i - N E_{P(v, h)}[v] \nabla_b = P(h| v^i) - N E_{P(v, h)}[h] \nabla_W = v^i P(h| v^i) - N E_{P(v^i, h)}[v^T h]

Correcciones

Usaremos como siempre el descenso de gradiente para corregir los casos, en este caso al ser un porblema de maximización no deberemos negar el gradiente.
a = a + lr \nabla_a b = b + lr \nabla_b W = W + lr \nabla_W

La gran problemática

Hemos deducido los gradientes de cada parámetro pero... M#@!¿a, volvemos a tener por allí la probabilidad conjunta de v y h que es extremadamente difícil de calcular... Como explicábamos al inicio de este mismo post, nos va a ser impracticable calcular la probabilidad conjunto y probablemente muramos antes obtener resultados.

La solución está en usar Gibbs sampling dentro del calculo del gradiente, vamos a aproximar P(v, h) a partir de las probabilidades condicionadas: P(v | h) y P(h | v). Para los que estéis familiarizados con la física probablemente lo habrías usado antes. Además vamos a hacer un engaño y decir que E_{P(v^i, h)}[x] \approx x \sim P(\hat{v}, \hat{h} | a, b, W) (la esperanza de todos los valores de x será una muestra de la estimación de la distribución de la probabilidad conjunta).
Como muchas veces sucede en este campo es mejor una aproximación rápida que una precisión brutal que tarde la vida, por varios motivos:

• Cuando tenemos un método iterativo, cuanto más rápido sea por iteración, más iteraciones realizaremos en el mismo marco de tiempo. Eso querrá decir que posiblemente estaremos más cerca de la solución optima con respecto otro método también iterativo donde su iteración sea mucho más costosa pero más exacta. Hay que buscar un compromiso.
• Trabajamos en entornos con muchísimo ruido, normalmente es más interesante priorizar tiempo que exactitud.

Gibbs sampling es un algoritmo de Markov chain Monte Carlo que permite iterativamente acercarse a la probabilidad conjunta real a costa de jugar al gato y al ratón con las probabilidades condicionadas. Permitiendo así ir acotando el resultado.

¿Que ganamos con esto? Aproximar {P(v^i, h)} lo suficiente para que el algoritmo sea capaz de aprender. A esto se le conoce como Contrastive Divergence.

¡Al código!

Para este ejemplo usaremos la base de datos MNIST junto a Tensorflow. Veréis como en ningún momento usamos Tensorflow para que nos realicé el proceso de optimización, sino que lo debemos hacer a mano al añadir el Gibbs sampling.

Para asegurarnos que el resultad es correcto, podemos usar P(v | h) para conocer la reconstrucción de una representación en concreto.

Hasta aquí este post sobre RBM, como podéis ver ha sido bastante completo y hasta hemos tirado de algoritmos estocásticos para hacer ciertas partes computables.

Referencias:

• https://en.wikipedia.org/wiki/Restricted_Boltzmann_machine
• http://deeplearning.net/tutorial/rbm.html
• Pro Deep Learning with Tensorflow [Libro]