Реализация функции активации softmax для нейронных сетей

Я использую функцию активации Softmaxна последнем слое нейронной сети. Но у меня проблемы с безопасной реализацией этой функции.

Наивной реализацией будет следующая:

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();

Это не очень хорошо работает для > 100 скрытых узлов, потому что во многих случаях y будет NaN(если y(f) > 709, exp (y(f)) вернет inf). Я придумал такую ​​версию:

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();

где safeExpопределяется как

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

Эта функция ограничивает ввод exp. В большинстве случаев это работает, но не во всех случаях, и мне не удалось выяснить, в каких случаях это не работает. Когда у меня есть 800 скрытых нейронов в предыдущем слое, это вообще не работает.

Однако, даже если бы это работало, я как-то «искажал» результат ИНС. Можете ли вы придумать какой-либо другой способ вычисления правильного решения? Существуют ли какие-либо библиотеки или приемы C++, которые я могу использовать для расчета точного результата этой ANN?

изменить:Решение, предложенное Итамаром Кацем, таково:

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();

И оно действительно математически то же самое. Однако на практике некоторые небольшие значения становятся равными 0 из-за точности с плавающей запятой. Интересно, почему никто никогда не записывает эти детали реализации в учебники.