log(1 - softmax(X))? Let $\vec X$ be a vector. The $\vec V = \mathrm{logsoftmax}(\vec{X})$ function is defined as:
$$v_i = \ln\left(\frac{e^{x_i}}{\sum_i e^{x_i}}\right)$$
This is provided in machine learning numerical packages for numerical stability.
Is there a numerically stable implementation of:
$$\ln\left(1 - \frac{e^{x_i}}{\sum_i e^{x_i}}\right)$$
provided in standard packages (e.g. PyTorch, or scipy, etc.)? What's a good way of computing this?
For example, if $e^{x_i}$ represents the (non-normalized) probaiblity of a label, then this is the log of the probability of an incorrect label.
 A: Usually these values are not computed alone: the entire collection of $v_i$ and $\log(1 - \exp(v_i))$ is needed.  That changes the analysis of computational effort.
To this end, let
$$\bar x = \log\left(\sum_{j} e^{x_j}\right) = x_k + \log\left(1 + \sum_{j\ne k} e^{x_j-x_k}\right)$$
for any index $k.$  The right hand expression shows how $\bar x - x_k$ can be computed in a numerically stable way when $k$ is the index of the largest argument, for then the argument of the logarithm is between $1$ and $n$ (the number of the $x_i$) and the sum can be accurately computed using exp (especially when the $x_j$ are ordered from smallest to largest in the sum).
The relation
$$\eqalign{
\log\left(1 - \frac{e^{x_i}}{\sum_{j} e^{x_j}}\right) &= \log\left(\frac{\sum_{j\ne i} e^{x_j}}{\sum_{j} e^{x_j}}\right) \\
&= \log\left(\sum_{j\ne i} e^{x_j}\right) - x_k + \log\left(\frac{e^{x_k}}{\sum_{j} e^{x_j}}\right) \\
&= \log\left(-e^{x_i} + \sum_{j} e^{x_j}\right) - x_k + v_k \\
&= \log\left(1 - e^{x_i - \bar x}\right) + (\bar x - x_k) + v_k
}$$
reduces the problem to finding that last logarithm, which is accomplished by applying a log1mexp function to the difference $x_i - \bar x.$
For $n$ arguments $x_i,$ $i=1,2,\ldots, n,$ the total effort to compute all $2n$ values is


*

*$n$ computations of the $v_i$ using logsoftmax.

*One computation of $\bar x - x_k$ (using $n-1$ exponentials and a logarithm).

*$n$ invocations of log1mexp.
A: One option is to use the numerically stable log-softmax implementation in combination with a numerically stable $\text{log1m_exp}(x):=\log(1-\exp(x))$ function.
I believe the following is pretty good for $\text{log1m_exp}$:


*

*if x > -0.693147 you use $\log(-\text{expm1}(x))$, 

*otherwise $\text{log1m}(\exp(x))$,


where $\text{log1m}(x):= \text{log1p}(-x)$. $\text{log1p}(x):=\log(1+x)$ is usually implemented (e.g. in numpy or base R, similarly the inverse function for $\text{log1p}$, i.e. $\text{expm1}(x) = \text{log1p}^{-1}(x)$).
A: Here is a possible way to do this, in Julia code (which I think is quite readable even if one does not know Julia):
function log1msoftmax(x::AbstractArray; dims=1)
  m = maximum(x; dims=dims)
  e = exp.(x .- m)
  s = sum(e; dims=dims)
  return log.((s .- e) ./ s)
end

This has the same complexity as a logsoftmax (it does one exp and one log call per entry).
