So when we usually linearize general relativity with respect to metric perturbations $g_{\mu\nu}\rightarrow g_{\mu\nu}+h_{\mu\nu}$, we compute the correction to the inverse of the metric to first order in $h$:$$g^{\mu\nu}\rightarrow g^{\mu\nu}-g^{\mu\rho}g^{\nu\tau}h_{\rho\tau}:=g^{\mu\nu}+h^{\mu\nu}$$ where we define $h^{\mu\nu}$ to be $h_{\mu\nu}$ with indexes lifted using the inverse of the background metric.

To get this result we ask that to the first order $$(g_{\mu\tau}+h_{\mu\tau})(g^{\tau\nu}+h^{\tau\nu})=\delta_{\mu}^{\nu}$$ Imposing that this holds exactly we get: $$(g_{\mu\tau}+h_{\mu\tau})h^{\tau\nu}=-h_{\mu\tau}g^{\tau\nu}$$ Inverting the first factor we have $$h^{\rho\nu}+h^{\rho\mu}h_{\mu\tau}g^{\tau\nu}=h^{\rho\mu}(\delta^{\nu}_{\mu}+h_{\mu\tau}g^{\tau\nu})=-g^{\rho\mu}g^{\nu\tau}h_{\mu\tau}$$ but I don't know how to solve this. I should invert $(\delta^{\nu}_{\mu}+h_{\mu\tau}g^{\tau\nu})$; is there a symbolic way to get to the result without using the explicit formula of the inversion of a matrix? (or equivalently: is perhaps the resulting expression simple?) Even better: is there any other (more or less physical) reasoning to get to the exact correction to the inverse metric?