That’s right. As long as the force (gradient of negative log density) you’re using is consistently clipped the same way, then you immediately get reversibility (i.e., detailed balance), and thus Metropolis-within-Gibbs for the Hamiltonian step is doing the right thing. Like regular leapfrog, it doesn’t follow the true gradient flow exactly, but that’s not required for the algorithm.
This is slightly changing the present topic, but if bisoftplus is a good idea, then shouldn’t a similar replacement for the logistic function also be a good idea? Suppose that we map [0, 1] \to \mathbf{R} using the function
y = \frac14\frac{x - \tfrac12}{x - x^2}.
We calculate the inverse from the quadratic formula:
x = \frac12\left(1 - \frac1{4y} \pm \sqrt{1 + \frac{1}{(4y)^2}}\right)
We use + when y > 0 and - when y < 0. There’s a numerical cancellation problem near y = 0, but perhaps that can be overcome. If so, then this shares with bisoftplus the favorable property that we approach the asymptotes only polynomially fast and not exponentially.
(Update: Come to think of it, we can eliminate the need to fiddle with the sign by factoring 1/(4y)^2 out of the square root; the resulting 1/(4y) has the requisite sign. But something more is necessary for numerical stability.)
Here’s a method for computing the inverse of my previous comment in a numerically accurate way. I don’t know whether my proposed function is statistically desirable, but at least it’s computationally feasible.
We use Newton’s method. Set
y = f(x) = \frac14 \frac{x - \tfrac12}{x(1 - x)}.
We fix some y_0 near zero and seek a root of f(x) - y_0 using Newton’s method. We compute
f'(x) = \frac14 \frac{x^2 - x + \tfrac12}{(x(1 - x))^2}.
We initialize x^{(0)} = y_0 + \tfrac12 and iterate with the usual equation
x \leftarrow x - \frac{f(x) - y_0}{f'(x)} = \frac{x(1-x)}{x^2 - x + \tfrac12}\big({-4x}(1 - x)y_0 + (x - \tfrac12)\big).
I tried it out, and it works fine as long as -\tfrac12 < y < \tfrac12. When \lvert y\rvert < 0.001 and we use double-precision arithmetic, we need only a single Newton step to reach our maximum precision.
I tried going up to second order and using Halley’s method. Some miraculous cancellation occurs and we get the update rule
x \leftarrow x - 4(x^2 - x + 1/2)\big({-4x}(1 - x)y_0 + (x - \tfrac12)\big).
So the update rule is strictly faster because it avoids a division. However, as is typical, this update rule doesn’t get really fast until we’re quite close, so when 0.0001 < \lvert y\rvert < 0.001 and we start at y + \tfrac12, this update rule needs two steps to reach our maximum precision; only when \lvert y\rvert < 0.0001 do we need just a single step.
Altogether, it looks to me like a good implementation could do:
- If \lvert y\rvert \ge 0.5, use the quadratic formula \tfrac12\big(1 + (4y)^{-1}(-1 + \sqrt{(4y)^2 + 1})\big).
- If 0.001 \le \lvert y\rvert < 0.5, initialize with the inverse from the quadratic formula and apply a single Newton step.
- If 0.0001 \le \lvert y\rvert < 0.001, initialize with x_0 = y + \tfrac12 and apply a single Newton step.
-
- If \lvert y\rvert < 0.0001, initialize with x_0 = y + \tfrac12 and apply a single Halley step.
But I’m not a numerical analyst by training, so it would be wise to run this by one and see what they think.
I don’t have a good sense of how this would change regressions, for example. Someone like @bgoodri or @avehtari probably would. Any cdf can be used here conceptually. Logistic regression uses the cdf of the standard logistic distribution, which is just the \textrm{logit}^{-1} function, aka the “logistic sigmoid.” Probit just uses the standard normal cdf.
The logistic function has an interpretation as the inverse of the log odds function. The log odds has a lot of nice properties beyond interpretability. For example, log odds are symmetric around 0 and thus have the nice property that \text{logit}^{-1}(-x) = 1 - \text{logit}(x). The new proposal also bakes in this symmetry with the splitting of inputs at 0.
Here’s a plot of the proposal vs. inverse logit (standard logistic cdf) and probit (standard normal cdf). Note that the x axis goes from -10 to 0 and the y axis is on the log scale to compare the tails.
Here’s one for x from -2 to 0 to get a better sense of how they start.
@kyle.hofmann You can certainly do things like that, and I think there are sometimes good reasons for doing so. If we apply the same asinh trick to a logistic regression, and use 1/2 (1 + tanh(asinh(x))) as link function in a regression, we also push predictions away from the extremes zero and one, just like we push predicted standard deviations away from zero and inf with the bisoftplus function. If my math is correct, this happens to be equivalent to a robit model with nu=2.
I see less of a reason to use this instead of a logit link in regressions by default though. In a a regression on the standard deviation with an exp link, we end up with a double exponential term like exp(exp(x)) in the density, so we still have an exp(x) term in the log-density, and that leads to problems (both computationally, but I think usually also model-wise). That’s not the case in a logistic regression.
I’m still thinking about this a little, but the bisoftplus function seems to have very nice properties specifically when talking about variances. There is something inherently additive about both variance and precision: if we add two independent random variables, the variances add up. And if we add measurements, the precision adds up. Now, let’s model a standard deviation as \sigma = \sigma_0 \text{bisoftplus}(\eta/2), where \sigma_0 is some estimated base-standard deviation and eta is a linear predictor without intercept. Then if \eta is close to zero, it works like an exp link up to a factor of 2. If eta is large, we get \sigma^2 = 2\sigma_0^2 + \eta^2 \sigma_0^2, so we can interpret eta as an additional source of variance, that we add the the base variance. And similarly for negative eta, only that eta acts like a measurement that makes the base variance smaller. In contrast, the exp link is all about multiplicative behaviour, which I don’t think is all that natural for variances.
I guess I need to write this down in a bit more detail, and add an example, but I think this really should be quite nice.
I agree with your reasoning for why bisoftplus might be desirable, but I am not so sure that that reasoning implies that my proposal is undesirable. What I hear you saying is that you can think of a natural statistical situation (regression on the standard deviation) where an exponential link creates an undesirable double exponential term in the density. Therefore a replacement link with non-exponential behavior is of real practical interest, and bisoftplus looks like it additionally has good additive behavior.
I also agree that, lacking a concrete situation where a logistic link creates a similarly problematic double exponential, there is less motivation for replacing it. For that reason, I can only refer to the link I proposed as an idea or a suggestion and not as a recommendation.
That said, here is a heuristic handwavy justification: Consider that for 0 < x \ll 1, the link I proposed has y \sim O(-1/x), and similarly for 0 \ll x < 1, it has y \sim O(1/(1-x)). Therefore the tails of the inverse link decay like O(1/y), which is even slower than the Cauchy distribution; you could call it “obese-tailed”. (I haven’t tried to check, but I wouldn’t be surprised if the inverse link defines a distribution with an undefined p’th moment for any p > 0.) (Update: I checked this. SymPy took a symbolic antiderivative for me and got the hypergeometric function {}_2F_1. Upon looking up its asymptotics, it was clear that none of the p’th moments were defined.) On a practical level, this means that an interval like (\epsilon, 1 - \epsilon) \subseteq [0, 1] corresponds to a very wide region of \mathbf{R}, something like (-\epsilon^{-1}, \epsilon^{-1}). Since the middle of [0, 1] is sent to such a wide interval, we can think of this link as keeping our attention focused on the middle of [0, 1] and away from the endpoints. For instance: Numbers as small as 10^{-100} don’t turn up in practice, so a good link function would make them correspond to points of \mathbf{R} that are so large we don’t expect them to turn up in practice. The fact that 10^{-100} \in [0, 1] corresponds to -1.25 \cdot 10^{99} under the obese-tailed link, to -230.26 under the logistic link, and to -21.27 under the probit link is a point in favor of the obese-tailed link and against the other two.
Of course, a successful application would be more important than a heuristic. To me, it sounds like bisoftplus has a promising application. If that works out, then perhaps this other link is worth some real effort.
I never said I that I think your link function is undesirable, only that I don’t think it would be a good default replacement of the logit link.
Other link functions for logistic regressions are used regularly, for instance the robit regression I mentioned. It is often framed as being robust to misclassifications, but you can just as well interpret it as a link function that give less prior weight on extreme probabilities close to 0 and 1. As I mentioned earlier, (1 + \tanh(\text{asinh}(x)) / 2 is equivalent to robit with \nu = 2, and similarly we can interpret your link function as a rational approximation to robit with \nu = 1. Which is also equivalent to a regression with cauchit link.
The cauchit link function is the inverse CDF function of a cauchy distribution, so g(p) = \tan(\pi(p - 1/2)) = -\cot(\pi p). We can find rational approximations of this with a bit of complex analysis (I don’t get to use that often, always nice when it comes up :-) ). This is a meromorphic function with poles at all integer values of p. We can write it using the Mittag-Leffler expansion – where each term corresponds to one pole – as
g(p) = -\frac{1}{\pi}\sum_{n=-\infty}^{\infty} \frac{1}{p-n}.
But we really care about the poles at 0 and 1, not all the others, so we just use those:
g(p) \approx \frac{1}{\pi}\left(\frac{1}{1 - p} - \frac{1}{p}\right) = \frac{1}{\pi}\frac{2p - 1}{p - p^2},
which is a rescaled version of your link function.
By construction from the Mittag-Leffler expansion those two have the same limiting behavior for p\to 0 and p \to 1. Around zero they differ a little bit:
My apologies; I seem to have misinterpreted your intent before. Also, thank you for the correction: I see now that I confused the rate at which the CDF approached 0 and 1 for the rate of decay of the PDF.
It is quite nice that the function I proposed has such a pretty relation to the Cauchy distribution! It makes me happy. So yes, what I am suggesting is nearly identical to using a Cauchy link. But a Cauchy link is certainly easier to use (since it’s already implemented) and better understood (since it’s a standard construction), and those are big advantages over my suggestion.
Thank you for the conversation. I look forward to hearing how bisoftplus performs in practice some day!