# Easiest way to have an approximation of binary entropy inverse

I call $H_2$ the function $H_2(x) = -x \ln(x)-(1-x) \ln(1-x)$

It is the binary entropy.

I call $g(x)=H_2((1+x)/2)$

This last function is bijective on $[0;1]$.

I would like to have an inverse, or at least an approximation of its inverse (because from what I understood, it is a transcendental equation to find its inverse).

I tried to use the function InverseFunction[g] but it didn't work (it doesn't compute, it returns $g^{-1}(0.5)$ and not its value if I apply it on $y=0.5$ for example).

How could I do it in the simplest way in mathematica ?

I like to use NDSolve to find inverses when symbolic methods won't work. The ODE corresponding to the inverse of a function g can be obtained by differentiating the defining equation:

eqn = g[ginv[h]] == h;
D[eqn, h]


g'[ginv[h]] ginv'[h] == 1

So, to find the inverse using NDSolve, we need the above equation and an initial point. For example:

ginv = NDSolveValue[
{ig'[h] == 1/g'[ig[h]], ig[g[1/2]] == 1/2},
ig,
{h, 0, 1},
WorkingPrecision->20
];


NDSolveValue::ndsz: At h == 0.6931471789228062290963080245570458567820., step size is effectively zero; singularity or stiff system suspected.

Plot[{g[x], ginv[x]}, {x, 0, 1}] I up the WorkingPrecision so that complex numbers don't pollute the result when the inverse approaches zero.

How about this? This samples the function g 100000 times an approximates the inverse of g by Interpolation:

H2 = x \[Function] -x Log[x] - (1 - x) Log[1 - x];
g = x \[Function] H2[(1 + x)/2];
xlist = Subdivide[0., 1., 100000];
ylist = Join[g[Most[xlist]], {0.}];
ginv = Interpolation[Transpose[{ylist, xlist}]];


Here is a plot of the approximate inverse function:

Show[
Plot[g[x], {x, xlist[], xlist[[-1]]}],
Plot[ginv[y], {y, ylist[], ylist[[-1]]}, PlotStyle -> ColorData],
PlotRange -> All, AspectRatio -> 1
] And this helps to estimate the error (it seems to be of order 10^-8):

Plot[{ginv[g[x]] - x, g[ginv[x]] - x}, {x, xlist[], xlist[[-1]]}, PlotRange -> All] Edit

With g defined as a pure function as above, you can also use ginv = InverseFunction[g]:

Plot[InverseFunction[g][y], {y, ylist[], ylist[[-1]]}] But plotting takes way longer, probably because InverseFunction involves a root finding algorithm for each function evaluation.

• Thank you for your answer. Indeed we can do trick but I wondered if there were like a "magic" function already implemented. But apparently not. I will use your trick (I'm a mathematica beginner so your example will be helpful) Apr 23 '18 at 14:21
• As I added in the post, InverseFunction does work for pure functions. Apr 23 '18 at 14:27
• Actually I don't understand how InverseFunction knew on which part he took the inverse function. Indeed the function is not bijective and mathematica chose one of the branches. But where has this choice been put in the code ? Apr 25 '18 at 9:50
• Honestly, I cannot tell what happens in the internals of InverseFunction... I guess it uses a root finding algorithm by bracketing the solution, but that's merely a shot into the dark. Apr 25 '18 at 9:55
• "Forced" is a strong word. The interpolation method may be more robust and also more efficient, so there is some reason to prefer it. See also Carl Woll's answer; there you should also be able to pick a branch for your choice. Apr 25 '18 at 10:00

This is the same approach as Henrik's answer (didn't notice the answer while I was working on it). Since my method to generate the dataset used for interpolation seems slightly more straightforward, I'll leave this here anyway:

Clear[h2, g]
h2[x_] := -x Log[x] - (1 - x) Log[1 - x]
g[x_] := h2[(1 + x)/2]


Generate an interpolation of the inverse by sampling g function:

ginv = Interpolation[
With[{pitch = 1/1000}, Table[{g[x], x}, {x, 0, 1 - pitch, pitch}]]
];


Plotting the results:

Show[
Plot[Legended[g[x], "g"], {x, 0, 1}, PlotRange -> All],
Plot[
Legended[ginv[x], "inverse"], Flatten@{x, ginv["Domain"]},
PlotStyle -> ColorData
],
Graphics[{Dashed, InfiniteLine[{{0, 0}, {1, 1}}]}],
AspectRatio -> 1
]
` 