Why $\delta R_{\mu\nu} = g^{\rho\sigma}\delta R_{\rho\mu\sigma\nu}$? Varying the Ricci tensor with respect to the metric $g^{\mu\nu}$, one would get
$
\delta R_{\mu\nu} = \delta(g^{\rho\sigma} R_{\rho\mu\sigma\nu}) = g^{\rho\sigma}\delta R_{\rho\mu\sigma\nu} + \delta g^{\rho\sigma} R_{\rho\mu\sigma\nu},
$
but in all my references I found that $\delta R_{\mu\nu} = g^{\rho\sigma}\delta R_{\rho\mu\sigma\nu}$ (including Wikipedia), which implies that the last term of the above equation is identically zero. Why is that so?
 A: Consider a metric $\widetilde{g}_{\mu\nu}=g_{\mu\nu}+\delta g_{\mu\nu}$ and its inverse metric $\widetilde{g}^{\mu\nu}=g^{\mu\nu}-\delta g^{\mu\nu}$.
(Fittingly I've just explained the negative sign here yesterday.) We have:
\begin{align*}
R_{\mu\nu}[\widetilde{g}]
&=\widetilde{g}^{\kappa\lambda}
R_{\kappa\mu\lambda\nu}[\widetilde{g}]
=(g^{\kappa\lambda}-\delta g^{\kappa\lambda})
(R_{\kappa\mu\lambda\nu}[g]+\delta R_{\kappa\mu\lambda\nu}) \\
&=R_{\mu\nu}[g]
+\underbrace{g^{\kappa\lambda}\delta R_{\kappa\mu\lambda\nu}
-\delta g^{\kappa\lambda}R_{\kappa\mu\lambda\nu}[g]}_{=\delta R_{\mu\nu}}+\mathcal{O}(\delta g^2).
\end{align*}
I therefore don't think this formula is correct in a stricter sense, somebody probably approximated $\widetilde{g}^{\kappa\lambda}\approx g^{\kappa\lambda}$ directly in the first step and therefore left out a term in first order. Because of this I'd also use a different formula for the pertubation of the curvature tensor. Taking the formula for the Christoffel symbols and putting the perturbed metric in, you get:
\begin{equation}
\Gamma_{\mu\nu}^\kappa[\widetilde{g}]
=\Gamma_{\mu\nu}^\kappa[g]
+\underbrace{\frac{1}{2}g^{\kappa\lambda}\left(
\nabla_\mu\delta g_{\lambda\nu}
+\nabla_\nu\delta g_{\lambda\mu}
-\nabla_\lambda\delta g_{\mu\nu}\right)}_{=\delta\Gamma_{\mu\nu}^\kappa}
+\mathcal{O}(\delta g^2).
\end{equation}
This calculation is not short, but pretty straight forward. You immediately see, that $\delta\Gamma_{\mu\nu}^\kappa$ is a tensor, as it is a sum of covariant differentiations of a tensor. Taking the formula for the Riemann curvature tensor and putting this equation in, you get:
\begin{equation}
R_{\lambda\mu\nu}^\rho[\widetilde{g}]
=R_{\lambda\mu\nu}^\rho[g]
+\underbrace{\partial_\mu\delta\Gamma_{\lambda\nu}^\rho
-\partial_\nu\delta\Gamma_{\lambda\mu}^\rho
+\delta\Gamma_{\mu\sigma}^\rho\Gamma_{\lambda\nu}^\sigma[g]
+\Gamma_{\mu\sigma}^\rho[g]\delta\Gamma_{\lambda\nu}^\sigma
-\delta\Gamma_{\nu\sigma}^\rho\Gamma_{\lambda\mu}^\sigma[g]
-\Gamma_{\nu\sigma}^\rho[g]\delta\Gamma_{\lambda\mu}^\sigma}_{=\delta R_{\lambda\mu\nu}^\rho}
+\mathcal{O}(\delta g^2).
\end{equation}
Notice, that you can shorten this using the covariant differentiation:
\begin{equation}
\delta R_{\lambda\mu\nu}^\rho
=\nabla_\mu\delta\Gamma_{\lambda\nu}^\rho
-\nabla_\nu\delta\Gamma_{\lambda\mu}^\rho.
\end{equation}
You again immediately see, that $\delta R_{\lambda\mu\nu}^\rho$ is a tensor, as it is the sum of covariant differentiations of a tensor. Contracting to $\delta R_{\lambda\nu}$ gives you the simple Palatini identity.
A: This is for whoever is interested in the solution. The second term in the equation of my original post is not zero as it is. But, from $g_{\mu\alpha}g^{\alpha\nu}=\delta_\mu^\nu$, one can prove that $\delta g_\rho^\sigma = - g_\alpha^\sigma g_\rho^\beta \delta g_\beta^\alpha$. Then the variation of the Ricci tensor gives
$
\delta R_{\mu\nu} = \delta(g_\rho^\sigma R^\rho_{\;\mu\sigma\nu})= g^\sigma_\rho\delta R^\rho_{\;\mu\sigma\nu} + \delta g_\rho^\sigma R^\rho_{\;\mu\sigma\nu},
$
where
$
\delta g_\rho^\sigma R^\rho_{\;\mu\sigma\nu} = - g_\alpha^\sigma g_\rho^\beta \delta g_\beta^\alpha R^\rho_{\;\mu\sigma\nu} = -\delta g_\beta^\alpha R^\beta_{\;\mu\alpha\nu},
$
which vanishes as $\rho,\sigma,\alpha,\beta$ are all dummy indices. Hence,
$
\delta R_{\mu\nu} = g^\sigma_\rho\delta R^\rho_{\;\mu\sigma\nu}.
$
A: The variation of the Ricci tensor
$$
\delta R_{\mu\nu}~=~\delta(g^{\rho\sigma} R_{\rho\mu\sigma\nu}) = g^{\rho\sigma}\delta R_{\rho\mu\sigma\nu} + \delta g^{\rho\sigma} R_{\rho\mu\sigma\nu}~=~g^{\rho\sigma}\delta R_{\rho\mu\sigma\nu}~+~\delta s\left(\frac{D}{ds} g^{\rho\sigma}\right)R_{\rho\mu\sigma\nu}.
$$
The last term is zero by the covariant constancy of the metric. This is seen with
$$
\frac{D}{ds} g^{\rho\sigma}~=~\frac{d}{ds} g^{\rho\sigma}~+~\left(\Gamma^\rho_{\alpha\beta}g^{\alpha\sigma}~+~\Gamma^\sigma_{\alpha\beta}g^{\alpha\rho}\right)\frac{dx^\beta}{ds}. 
$$
This is why the $\delta g^{\rho\sigma} R_{\rho\mu\sigma\nu}~=~0$
