Implementación de una función de activación softmax para redes neuronales

Question

Estoy usando una función de activación Softmax en la última capa de una red neuronal. Pero tengo problemas con una implementación segura de esta función.

Una implementación ingenua sería éste:

Vector y = mlp(x); // output of the neural network without softmax activation function 
for(int f = 0; f < y.rows(); f++) 
    y(f) = exp(y(f)); 
y /= y.sum(); 

Esto no funciona muy bien para> 100 nodos ocultos debido a la y será NaN en muchos casos (si y (f)> 709, exp (y (f)) devolverá inf). Se me ocurrió con esta versión:

Vector y = mlp(x); // output of the neural network without softmax activation function 
for(int f = 0; f < y.rows(); f++) 
    y(f) = safeExp(y(f), y.rows()); 
y /= y.sum(); 

donde safeExp se define como

double safeExp(double x, int div) 
{ 
    static const double maxX = std::log(std::numeric_limits<double>::max()); 
    const double max = maxX/(double) div; 
    if(x > max) 
    x = max; 
    return std::exp(x); 
} 

Esta función limita la entrada de exp. En la mayoría de los casos, esto funciona, pero no en todos los casos, y realmente no pude averiguar en qué casos no funciona. Cuando tengo 800 neuronas ocultas en la capa anterior, no funciona en absoluto.

Sin embargo, incluso si esto funcionó, de alguna manera "distorsionar" el resultado de la ANN. ¿Puedes pensar en otra forma de calcular la solución correcta? ¿Hay alguna biblioteca C++ o trucos que pueda usar para calcular la producción exacta de esta ANN?

edición: La solución proporcionada por Itamar Katz es:

Vector y = mlp(x); // output of the neural network without softmax activation function 
double ymax = maximal component of y 
for(int f = 0; f < y.rows(); f++) 
    y(f) = exp(y(f) - ymax); 
y /= y.sum(); 

Y realmente es matemáticamente la misma. En la práctica, sin embargo, algunos valores pequeños se vuelven 0 debido a la precisión del punto flotante. Me pregunto por qué nadie escribe estos detalles de implementación en los libros de texto.

Answer 1

Primero ve a la escala de registro, es decir, calcula log(y) en lugar de y. El registro del numerador es trivial. Para calcular el registro del denominador, puede usar el siguiente 'truco': http://lingpipe-blog.com/2009/06/25/log-sum-of-exponentials/

Answer 2

Sé que ya está respondida, pero la publicaré paso a paso de todos modos.

poner en el registro:

zj = wj . x + bj 
oj = exp(zj)/sum_i{ exp(zi) } 
log oj = zj - log sum_i{ exp(zi) } 

Sea M el MAX_I {zi} utilizar el truco de registro de suma-exp:

log oj = zj - log {sum_i { exp(zi + m - m)}} 
    = zj - log {sum_i { exp(m) exp(zi - m) }}, 
    = zj - log {exp(m) sum_i {exp(zi - m)}} 
    = zj - m - log {sum_i { exp(zi - m)}} 

el término exp (zi-m) pueden sufrir desbordamiento si m es mucho mayor que otros z_i, pero está bien, ya que esto significa que z_i es irrelevante en la salida de softmax después de la normalización. los resultados finales son:

oj = exp (zj - m - log{sum_i{exp(zi-m)}})