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 "by-hand" solution using ReplaceRepeated agree.