This sum could be really large, and reach infinity, or at least be very imprecise, and I don't think the optimizations to mitigate that exist in Theano.

I think it would be better not to try and craft a `y` so that `softmax(y)` has the desired min_val, but instead to compute softmax(x) (or plushmax(x) with the original epsilon), and then scale that result between `min_val` and `min_val*(1-n_classes)`.
The formula for that would be `softmax(y) * (1 - n_classes * min_val) + min_val`, with `y` computed with the _original_ `eps`.

I'm not sure if it is really important not to do that scaling if `eps` is already large enough for `softmax(y).min()` to be above `min_val`, but if it is, we can always do `m = max(min_val, (min_val - softmax(y).min()) / (1 - n_classes * softmax(y).min()))` and use `softmax(y) * (1 - n_classes * m) + m`.