5
$\begingroup$

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 ?

$\endgroup$
7
$\begingroup$

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.69314717892280622909630802455704585678`20., step size is effectively zero; singularity or stiff system suspected.

Plot[{g[x], ginv[x]}, {x, 0, 1}]

enter image description here

I up the WorkingPrecision so that complex numbers don't pollute the result when the inverse approaches zero.

$\endgroup$
4
$\begingroup$

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[[1]], xlist[[-1]]}],
 Plot[ginv[y], {y, ylist[[1]], ylist[[-1]]}, PlotStyle -> ColorData[97][2]],
 PlotRange -> All, AspectRatio -> 1
 ]

enter image description here

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[[1]], xlist[[-1]]}, PlotRange -> All]

enter image description here

Edit

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

Plot[InverseFunction[g][y], {y, ylist[[1]], ylist[[-1]]}]

enter image description here

But plotting takes way longer, probably because InverseFunction involves a root finding algorithm for each function evaluation.

$\endgroup$
  • $\begingroup$ 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) $\endgroup$ – StarBucK Apr 23 '18 at 14:21
  • 2
    $\begingroup$ As I added in the post, InverseFunction does work for pure functions. $\endgroup$ – Henrik Schumacher Apr 23 '18 at 14:27
  • $\begingroup$ 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 ? $\endgroup$ – StarBucK Apr 25 '18 at 9:50
  • $\begingroup$ 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. $\endgroup$ – Henrik Schumacher Apr 25 '18 at 9:55
  • 1
    $\begingroup$ "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. $\endgroup$ – Henrik Schumacher Apr 25 '18 at 10:00
3
$\begingroup$

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:

Define your functions:

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[97][2]
 ],
 Graphics[{Dashed, InfiniteLine[{{0, 0}, {1, 1}}]}],
 AspectRatio -> 1
]

plot of g and its inverse

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.