I have found an answer using Solve. It's not a single function call as I'd hoped, but it could probably be made into a function. Perhaps other users will have improvements to this answer.
Consider H[L_, c] := ... as defined in the question. Define a series representation for the first-argument inverse, L[H,c]:
L[H, c] = L0[H] + 1/c^2 L2[H] + 1/c^4 L4[H]
Here only the even terms are required as can be seen in the definition of H.
Now we can setup the equation H == H[L[H, c], c] and ask Mathematica to Solve for the coefficients at each order in 1/c. For convenience, to avoid imaginary numbers in the solution, I have defined a dummy variable $\varepsilon = - 2 H$.
eqs = Map[# == 0 &, CoefficientList[(-ɛ/2) - H[L[H, c], c] + O[c, Infinity]^6, 1/c]];
Print /@ Select[eqs, ! MatchQ[#, True] &];
sol = Last@Solve[eqs, {L0[H], L2[H], L4[H]}] /. ɛ -> -2 H
Here, Last was used to take the positive solution, as discussed in the question. The final result is:
L[H, c] /. sol
It can be verified that this solution and the t"by-hand" solution using ReplaceRepeated agree.
Edit: As pointed out by @yarchik, SolveAlways is designed precisely for setting polynomial equations to zero and solving for coefficients. While I couldn't get SolveAlways to solve for L0, L2, and L4 directly, an equivalent formulation can be achieved using an Eliminate one-liner:
Solve[Eliminate[(-ɛ/2) == H[L[H, c], c] + O[c, Infinity]^6, c], {L0[H], L2[H], L4[H]}] // Simplify // ReplaceAll[ ɛ -> -2 H]