How to tell Mathematica to get the real-only inverse function?

Reduce[expr, x, Reals] will be your friend here, but it can take a bit of work to parse its result.

Here's a solution that should work for any expression, not just polynomials (at least for the small set of examples I tried).

RealInverse[a_. x_^q_Integer?Positive + b_., x_] /; FreeQ[{a, b}, x] := 
  Surd[(x-b)/a, q]

RealInverse[expr_, x_] := 
  Module[{y, red},
    red = Reduce[x == (expr /. x -> y), y, Reals];
    ToPW[PWTerm[#, y]& /@ OrList[BooleanConvert[red, "DNF"]], x]
  ]

OrList[HoldPattern[Or][args__]] := {args}
OrList[e_] := {e}

PWTerm[And[cond1___, y_ == root_, cond2___], y_] /; FreeQ[{root, And[cond1, cond2]}, y] := {root, And[cond1, cond2]}
PWTerm[y_ == root_, y_] /; FreeQ[root, y] := {root, True}
PWTerm[___] = $Failed;

ToPW[lis_?MatrixQ, x_] /; Length[First[lis]] == 2 := 
  With[{condlist = Last[Transpose[lis]]},
    Piecewise[lis, Undefined] /; DisjointConditionsQ[condlist, x]
]
ToPW[___] = $Failed;

DisjointConditionsQ[{_}, _] = True;
DisjointConditionsQ[cond_List, x_]:=
  Reduce[And @@ Table[DisjointCondition[cond, i], {i, Length[cond]}], x]

DisjointCondition[cond_, i_] := 
  If[TrueQ[cond[[i]]],
    True,
    Implies[cond[[i]], Not[And @@ Delete[cond, i]]]
]

Now some tests:

RealInverse[x^3, x]
Surd[x, 3]
RealInverse[x^7 + x^4 + 3 x - 1, x]
Root[-1 - x + 3 #1 + #1^4 + #1^7 &, 1]
(* Not invertible over the real line *)
RealInverse[x^3 - 3 x + 1, x]
$Failed
RealInverse[x^(2/3) + Sqrt[x], x]
ConditionalExpression[
  Root[x^6 - 3 x^4 #1 + (3 x^2 - 2 x^3) #1^2 + (-1 - 6 x) #1^3 + #1^4 &, 1], 
  x >= 0
]
RealInverse[Sign[x] (Abs[x]^(2/3) + Sqrt[Abs[x]]), x]

enter image description here