Perturbation metric problem

Note that: $$ h^{\mu \nu} = \eta^{\mu \rho}\eta^{\nu \lambda} h_{\rho \lambda} $$ Therefore, up to first order, we have: \begin{equation} \begin{aligned} g^{\mu \nu}g_{\nu \sigma} & = (\eta^{\mu \nu} - h^{\mu \nu})(\eta_{\nu \sigma} + h_{\nu \sigma}) \\& =\eta^{\mu \nu}\eta_{\nu \sigma} + \eta^{\mu \nu}h_{\nu \sigma} - \eta_{\nu \sigma} h^{\mu \nu} + \mathcal{O}(h^2) \\& = \delta^\mu_\sigma + \eta^{\mu \nu}h_{\nu \sigma} - \eta_{\nu \sigma} \eta^{\mu \rho}\eta^{\nu \lambda} h_{\rho \lambda} + \mathcal{O}(h^2) \\& = \delta^\mu_\sigma + \eta^{\mu \nu}h_{\nu \sigma} - \delta_\sigma^\lambda \eta^{\mu \rho} h_{\rho \lambda} + \mathcal{O}(h^2) \\& = \delta^\mu_\sigma + \eta^{\mu \nu}h_{\nu \sigma} - \eta^{\mu \rho} h_{\rho \sigma} + \mathcal{O}(h^2) \\& = \delta^\mu_\sigma + \eta^{\mu \nu}h_{\nu \sigma} - \eta^{\mu \nu} h_{\nu \sigma} + \mathcal{O}(h^2) \\& = \delta^\mu_\sigma + \mathcal{O}(h^2) \end{aligned} \end{equation} which is what we require.

Edit in response to comments. The inverse of $g_{\nu \sigma}$, which is denoted by $g^{\mu \nu}$, is defined to satisfy the following equation: $$ g^{\mu \nu}g_{\nu \sigma} = \delta^\mu_\sigma $$ In other words, we need to find an expression for $g^{\mu \nu}$ such that the following equations is satisfied: $$ g^{\mu \nu}(\eta_{\nu \sigma} + h_{\nu \sigma}) = \delta^\mu_\sigma $$ As I have shown above, the function that obeys the above equation up to first order in $h_{\mu \nu}$ is: \begin{equation} g^{\mu \nu} = \eta^{\mu \nu} - h^{\mu \nu} \end{equation} That is really all that is going on.

If you have $g_{\mu\nu} = \eta_{\mu\nu} + h_{\mu\nu}$, with $|h_{\mu\nu}|\ll 1$ a perturbation, $\eta_{\mu\nu} = \text{diag}(-1,1,1,1)$, then you can simply perform a Taylor expansion to obtain the inverse, $$g^{\mu\nu} = \eta^{\mu\nu} - h^{\mu\nu}+ \mathcal{O}(h^2),$$ where the indices of $h^{\mu\nu}$ are raised by $\eta^{\mu\nu}$. This is easy to see if you take the diagonal elements first, the off-diagonal ones are of second order, negligible. You don't need to solve a matrix equation for this, because the inverse is linearly dependent on $\eta_{\mu\nu}=\eta^{\mu\nu}$, $h_{\mu\nu}$ and their products, hence it must be of the form $g^{\mu\nu} = a\eta^{\mu\nu} + b h^{\mu\nu} + \ldots$, where the ellipses imply linear dependencies of higher order and the coefficients are constants. Also, note that $h_{\mu\nu} \eta^{\mu\mu} \eta^{\nu\nu} = h^{\mu\nu}$, no summation implied, namely $h^{\mu\nu} \propto h_{\mu\nu}$.

Also, a straightforward way to calculate $\delta g^{\mu\nu}$ is by use of $\delta g^{\mu\nu} = - g^{\mu\rho} g^{\nu\sigma} \delta g_{\rho\sigma}$ with $\delta g_{\mu\nu} = h_{\mu\nu}$.