As you said, I suppose that the function$$h(x) = -x \log_{2}(x) -(1-x) \log_{2}(1-x)\tag1$$ is not invertible.
However, if you are not too concerned by the values close to the boundaries, $h(x)$ can be approximated by a $[4,4]$ Padé approximant built at $x=\frac 12$. This would give
$$h(x) \sim \frac{1- \left(\frac{124}{49}+\frac{2}{\log
(2)}\right)\left(x-\frac{1}{2}\right)^2+ \left(\frac{152}{245}+\frac{548}{147 \log
(2)}\right)\left(x-\frac{1}{2}\right)^4 } {1-\frac{124}{49} \left(x-\frac{1}{2}\right)^2+\frac{152}{245}
\left(x-\frac{1}{2}\right)^4 }\tag2$$ which reduces to a quadratic equation in $\left(x-\frac{1}{2}\right)^2$.
For sure, this could be improved using the $[6,6]$ Padé approximant which would reduces to a cubic equation in $\left(x-\frac{1}{2}\right)^2$ which would be less pleasant.
In order to check the quality of the approximation, give $x$ a value in $(1)$ to get $h(x)$ and for this value, solve $(2)$ to get the solution.
$$\left(
\begin{array}{cc}
x_{\text{given}} & x_{\text{recomputed}} \\
0.05 & 0.048894 \\
0.10 & 0.099763 \\
0.15 & 0.149949 \\
0.20 & 0.199990 \\
0.25 & 0.249998 \\
0.30 & 0.300000 \\
0.35 & 0.350000 \\
0.40 & 0.400000 \\
0.45 & 0.450000 \\
0.50 & 0.500000 \\
0.55 & 0.550000 \\
0.60 & 0.600000 \\
0.65 & 0.650000 \\
0.70 & 0.700000 \\
0.75 & 0.750002 \\
0.80 & 0.800010 \\
0.85 & 0.850051 \\
0.90 & 0.900237 \\
0.95 & 0.951106
\end{array}
\right)$$
Edit
The approximation can be improved using again
$$h(x) \sim f(x)=\frac{1+ a\left(x-\frac{1}{2}\right)^2+ b\left(x-\frac{1}{2}\right)^4 } {1+c \left(x-\frac{1}{2}\right)^2+d \left(x-\frac{1}{2}\right)^4 }$$
In order to respect the boundary conditions, we need to fix $\color{red}{b=-4 (a+4)}$. Now, considering the norm
$$\Phi(a,c,d)=\int_0^1 \big( h(x)-f(x) \big) ^2\,dx $$ its numerical minimization with respect to $a$, $c$ and $d$ (this is equivalent to a nonlinear regression for an infinte number of data points) leads to
$$\{a=-7.14483221448716,\,\, c=-4.28551864360839,\,\,d=2.63456254971723\}$$ At this point $\Phi(a,c,d)=3.44\times 10^{-7}$ while using the Padé approximant the norm would be $3.94\times 10^{-5}$ (that is to say $114$ times larger).
Update
If we take the problem for any base $a$,we have
$$h_a(x) = -x \log_{a}(x) -(1-x) \log_{a}(1-x)$$ that is to say
$$h_a(x)\log_e(a)=-x \log_{e}(x) -(1-x) \log_{e}(1-x)$$ and the $[4,4]$ Padé approximant of the rhs is
$$\frac{\log_e(2) + a\left(x-\frac{1}{2}\right)^2+ b\left(x-\frac{1}{2}\right)^4 }{{1-\frac{124}{49} \left(x-\frac{1}{2}\right)^2+\frac{152}{245}
\left(x-\frac{1}{2}\right)^4 } } $$ where
$$a=-\left(2+\frac{124 }{49}\log_e (2)\right)\qquad \text{and} \qquad b=\frac{548}{147}+\frac{152 }{245}\log_e (2)$$ and the norm is equal to $1.89\times 10^{-5}$.
Update
Answering this old question, I proposed two rather good approximations of the inverse
$$x_0=\frac{1}{2} \left(1-\sqrt{1-h^{4/3}(x)}\right)$$
$$x_1=\frac{\log \left(1+\sqrt{1-h^{4/3}(x)}\right)+(y-1) \log (2)}{2 \tanh
^{-1}\left(\sqrt{1-h^{4/3}(x)}\right)}$$
If $k=b/a$ is not very small or very big then I'd parametrize in this way:
$$\begin{align}
p_1 &= p^{1+x} \\
p_2 &= p^{1-x} \\
x &= \frac{\log(p_1/p_2)}{\log(p_1 p_2)} = \frac{b-a
}{b+a} =\frac{k-1}{k+1} \\
p &= \sqrt{p_1 p_2}=\exp\left(- \frac{a+b}{2ab}\right)=
\exp\left(- \frac{1}{a} \frac{1}{1+x}\right)
\end{align}
$$
Notice that $p$ is somewhere between $p_1$ and $p_2$, and if (say) $\frac{1}{9}<k<9$ then $-0.8<x<0.8$.
Doing a Taylor expasion of $f$ around $x=0$, we get the approximation
$$ f \approx -\log{\left( 1-p\right) }+\frac{p \left( 2 {{p}}-1\right) {{\log^{2}{(p)}}} }{6 (1-p)^2}{{x}^{2}} $$
(I'm using natural logarithms here, entropy in nats - to get it in bits just divide everything by $\log(2)$)
The approximation seems quite good, assuming $|x|$ is not very near $1$ and $p$ is not too small. Here's a graph for three values of $p$
and $x\in [-0,8,0.8]$. The dashed lines correspond to the approximation.
If you need more precision you can add the next term, which is
$$\frac{p\, \left(7-28p +47 p^2 - 8 p^3 \right) {{\log^4{(p)}}}}{360 {{\left(1-p\right) }^{4}}} x^4$$
Best Answer
There is a more systematic treatment of this question which I can refer to here. So this "replaces" in a way the precise bound in the question, I hope this is fine with OP.
One can ask the following question: suppose we want to establish an upper bound $$H_2(x) \leq c \cdot (x(1-x))^{q} \quad ,$$ with some constant $c$, what is the largest $q^*$ where this is possible?
You have already noticed that $q= 0.5$ is possible. You tried to improve that and found that $q= 0.75$ is "almost" possible. To make it strict, you (and Jack D'Aurizio) needed an extra additive constant.
The following results are taken from the paper "Bounds for entropy and divergence for distributions over a two-element set" by Topsøe, Flemming. In: JIPAM. Journal of Inequalities in Pure & Applied Mathematics. 2 (2): Paper No. 25 (2001). Link to pdf. Note when you read the paper that Topsøe uses natural logarithms so a little transformation is necessary.
From Topsøe's paper, the answer to the above question is $q^* = 1 / \ln(4) \simeq 0.721$ (which is close to your $0.75$). In this case, $c^* = 4^{q^*} = e \simeq 2.718 $, which can be compared to the constant in your proposition, $c = 4^{3/4} \simeq 2.82$, again very close.
Since $c^* = 4^{q^*}$ this is also a very good reason to include a factor $4$ into the RHS of the proposition, which you and Topsøe have done. Another reason is establishing the equality at $x^* = 0.5$, where for all $q$ holds: $1 = H_2(x^*) = (4 \, x^*(1-x^*))^{q} $.
Again, I modified your question for the reason of improving the somewhat arbitrary constant + exponent proposition and to make the treatment more concise, hope that is fine with you.