The speed of gravitational waves in general relativity

What is the speed of a wave in a non-linear theory? Answering before considering your question is important, because that answer will tell you where to look for your answer.

A useful notion is that of domain of dependence (see for example a decent book on GR for a detailed discussion, e.g., Wald or Hawking & Ellis). If $S\subseteq \Sigma$ is a subset of a Cauchy surface, then the domain of dependence $D(S)$ is the region of the spacetime where the solution is completely determined by the initial data on $S$, irrespective of what initial data is specified on the complement $\Sigma \setminus S$. Thus, the "slope" of the boundaries of $D(S)$ when represented in a spacetime diagram can be interpreted as the rate at which the influence of the initial data from $\Sigma \setminus S$ is encroaching on the spacetime region where the solution is determined by the initial data on $S$ alone. In other words, the "slope" of the boundary of $D(S)$ determines the speed of the propagation of disturbances in the solutions of the PDE you are solving. This speed also agrees with the speed of traveling wave solutions in the case of a linear PDE with constant coefficients.

Note that this "speed" of propagation could be independent of the underlying solution itself (the case in electromagnetism) or may actually depend on the solution (the case in GR). It happens to be a fact of life that for both electromagnetism and GR that the boundary of $D(S)$ is always traced out by the null rays of the spacetime metric (which is dynamical in the latter case, but need not be in the former). This fact then implies that the speed of electromagnetic and gravitational waves is the same, independent of any linearization.

The above result on the boundary of the domain of dependence is standard in the hyperbolic PDE literature and can be found, for instance, in the books of John (Partial differential equations), Lax (Hyperbolic partial differential equations) and Hoermander (Lectures on nonlinear hyperbolic differential equations). Briefly, the boundary is always "characteristic", where for PDEs of Lorentzian wave equation type the "characteristics" coincide with the null directions of the metric.