A ridiculously fast way to find fixed points of the logistic function? For $c, d > 0$, suppose I want to solve for $t$ in
\begin{align*}
t = \sigma(c - dt)
\end{align*}
where $\sigma(x) = (1 + \exp(-x))^{-1}$ is the logistic sigmoid function, also known as the expit function in statistics. I want an incredibly fast algorithm to find this, something that delivers $10^{-6}$ precision with one iteration or so.
Background and Motivation
In statistics & machine learning, there are many kinds of so-called (inverse) link functions which map a linear combination of covariates/features into a subset of $\mathbb{R}$. For example $\exp: \mathbb{R}\rightarrow \mathbb{R}_+$ and $\sigma: \mathbb{R} \rightarrow (0, 1)$, common inverse link functions for positive data and binary probabilities, respectively. I'm working on a class of numerical methods which, at some step in the process, need to find the fixed point of the inverse link function extremely quickly. Let's first look at the easier problem
\begin{align*}
t = \exp(c - dt) 
\end{align*}
Replacing $z = dt$,
\begin{align*}
ze^{z} = de^{c}
\end{align*}
This functional equation is what defines the Lambert $W$ function, and therefore $z = W(de^c) \implies t = W(de^c)/d$. The Lambert $W$ function appears quite often in physics and complex analysis, and therefore has been studied very well and many methods have been implemented for this. (See, for example, apparently the fastest implementation for the Lambert $W$ function by Darko Veberic.)
I tried to mimic the derivation of the above implementation by Veberic, but encountered the following problems, mostly simply due to the fact that $\sigma$ does not have as nice of properties as $\exp$:


*

*To derive an initial approximation, Veberic performed a series expansion around $W^{-1}(z) = z e^z$ at its extrema, and then inverted. Unfortunately, such extrema do not exist with $\sigma$.

*Unlike the $\exp$, I cannot isolate my unknown variable to one side like $z e^z = de^c$.


Any help on this would be greatly appreciated!

As in the comments I write $f(t)=\sigma(c-dt)$.
You can get a pretty fast method by just doing one Newton step starting from a relatively good guess. Since we know any fixed point must be in $\mathrm{range}(f)=(0,1)$, it is reasonable to start with a fixed point of some simple global approximation of $f$ on $(0,1)$. One is $g(t)=f(0)+(f(1)-f(0))t$. The fixed point of $g$ is $\frac{f(0)}{1-f(1)+f(0)}$. You can use that as an initial guess for Newton or Halley iteration and it should work pretty well because $f$ is reasonably nice, particularly if $d$ is not much larger than $c$.
Better guesses exist, using the information given about $f$ in more detail.
