Sigmoid's stability

Analytically, the logarithm of the sigmoid is always defined and finite, because the sigmoid returns values restricted to the open interval (0, 1), rather than using the entire closed interval of valid probabilities [0, 1]. In software implementations, to avoid numerical problems, it is best to write the negative log-likelihood as a function of z, rather than as a function of ˆy = σ( z). If the sigmoid function underflows to zero, then taking the logarithm of ˆy yields negative infinity.

Questions:

1. Sigmoid function's range is closed interval [0,1]? https://upload.wikimedia.org/wikipedia/commons/thumb/8/88/Logistic-curve.svg/1200px-Logistic-curve.svg.png Why does it mention as open?
2. How does write function of z rather y affect it from not reaching 0? What are those forms?

1. The sigmoid range $(0,1)$ is technically open, because no input value maps to $0$ or $1$. You can get arbitrarily close to $0$ or $1$ but never equal them.

2. Taking an exponent followed by a log can cause overflow issues during computation, because exponents grow quickly. Sigmoid also gets close to zero quickly and this can underflow (potentially rounding to 0). However, the end result of the combined function is likely to be something within normal calculation range, because of the log (this is similar to - but even more extreme - multiplying and dividing by a very large number).

There is no analytical form where $E = -log(\sigma(z))$ can be re-written as a simple function of z. In theory, it can be calculated numerically by expanding terms carefully (beyond me, so won't show here), and some libraries may include this kind of expansion. In practice, to avoid numerical instability, many machine learning libraries will simply cap the value like this: $E = -log( max(\epsilon, \sigma(z)) )$ with $\epsilon$ a small number, maybe $10^{-15}$
A more common issue where the calculations are done accurately in neural networks is when considering the gradient of a loss function. If you use a sigmoid output layer alongside a binary cross-entropy cost function $E = -(y log(\sigma(z)) + (1-y)log(1-\sigma(z)))$, then some of the terms cancel out, and then the gradient contribution from training is trivially $\frac{\partial{E}}{\partial z} = y - \sigma(z)$ - there is no need in that case to calculate any log values despite them being in the loss function.
1. You are correct, the sigmoid's range is $(0,1)$. The quote you shared is just saying, trivially, that anything in $[0,1]$ is a valid probability.
2. The sigmoid function gets very small, very quickly (because $e^{-z}$ grows exponentially as $z \to - \infty$ by definition). So the computer will round $\sigma(z)$ to zero for values of $z$ that aren't even that large in magnitude. This means you shouldn't use the output $y$ of the sigmoid in the log-likelihood; you should just build the sigmoid into the likelihood in a way that avoids numerical issues. Then, modestly large inputs $z$ won't give you numerical problems.