How to calculate $\sigma^{-1} (\frac{1}{n} \sum_{i=1}^n \sigma(x_i))$ in a numerically stable way? Suppose I have $n$ logits $x_1, \dots, x_n$, where $n$ is not too large. They are real numbers in $(-\infty, +\infty)$ and correspond to probabilities in $(0, 1)$ via the formula $p_i = \sigma(x_i)$, where $\sigma$ is the sigmoid function. I want to calculate the logit $\hat x = \sigma^{-1} (\frac{1}{n} \sum_{i=1}^n p_i) = \sigma^{-1} (\frac{1}{n} \sum_{i=1}^n \sigma(x_i))$ in a numerically stable way. How do I do this? Ideally, I want a code snippet for pytorch.

Well, I suppose that $\sigma$ is just the sigmoidal, which Pytorch has for sure an implementation, so you just need to find the inverse of the sigmoidal:
$$
\left(\frac{1}{1+e^{-x}}\right)^{-1} = \log\left(\frac{x}{1-x}\right)
$$
Then you can do whatever you want with this to make it more stable, for example:
$$
\log\left(\epsilon_1 + \frac{x}{\epsilon_2 + 1-x}\right)
$$
to avoid division by zero or the log of zero
