Since nobody has done so yet, I will address your first question. I also could not find a reference for [a proof of] this result anywhere, so if anyone knows a reference please let us know.
The most general test that this $F$-test can handle is $H_0 : C \beta = \psi$ for some $q \times p$ matrix $C$ and $q$-vector $\psi$. This allows you to test hypotheses like $H_0 : \beta_1 + \beta_2 = \beta_3 + 4$.
However, it seems you are focusing on testing hypotheses like $H_0 : \beta_2 = \beta_4 = \beta_5 = 0$, which is a special case with $\psi=0$ and $C$ being a matrix with one $1$ in each row, and all other entries being $0$. This allows you to more concretely view the smaller model as obtained by simply dropping some columns of your design matrix (i.e. going from $X_u$ to $X_r$), but in the end the result you are seeking is in terms of an abstract $C$ anyway.
Since it happens to be true that the formula $(C\hat{\beta})' (C (X'X)^{-1} C') (C \hat{\beta})$ works for arbitrary $C$ and $\psi=0$, I will prove it in that level of generality. Then you can consider your situation as a special case, as described in the previous paragraph.
If $\psi \ne 0$, the formula needs to be modified to $(C\hat{\beta} - \psi)' (C (X'X)^{-1} C') (C \hat{\beta} - \psi)$, which I also prove at the end of this post.
First I consider the case $\psi=0$. I will try to keep some of your notation. Let $V_u = \text{colspace}(X) = \{X\beta : \beta \in \mathbb{R}^p\}$. Let $V_r := \{X\beta : C\beta=0\}$. (This would be the column space of your $X_r$ in your special case.)
Let $P_u$ and $P_r$ be the projections on these two subspaces. As you noted, $P_u y$ and $P_r y$ are the predictions under the full model and the null model respectively. Moreover, you can show $\|(P_u - P_r) y\|^2$ is the difference in the sum of squares of residuals.
Let $V_l$ be the orthogonal complement of $V_r$ when viewed as a subspace of $V_u$. (In your special case, $V_l$ would be the span of the columns of the removed columns of $X_u$.) Then $V_r \oplus V_l = V_u$, and moreover, In particular, if $P_l$ is the projection onto $V_l$, then $P_u = P_r + P_l$.
Thus, the difference in the sum of squares of residuals is $$\|P_l y\|^2.$$ If $\tilde{X}$ is a matrix whose columns span $V_l$, then $P_l = \tilde{X} (\tilde{X}'\tilde{X})^{-1} \tilde{X}'$ and thus $$\|P_l y\|^2 = y'\tilde{X} (\tilde{X}'\tilde{X})^{-1} \tilde{X}' y.$$
In view of your attempt at the bottom of your post, all we have to do is show that choosing $\tilde{X} := X(X'X)^{-1} C'$ works, i.e., that $V_l$ is the span of this matrix's columns. Then that will conclude the proof.
- It is clear that $\text{colspace}(\tilde{X}) \subseteq \text{colspace}(X)=V_u$.
- Moreover, if $v \in V_r$ then it is of the form $v=X\beta$ with $C\beta=0$, and thus $v' \tilde{X} = \beta' X' X (X'X)^{-1} C' = (C \beta)' = 0$, which shows $\text{colspace}(\tilde{X})$ is in the orthogonal complement of $V_r$, i.e. $\text{colspace}(\tilde{X}) \subseteq V_l$.
- Finally, suppose $X\beta \in V_l$. Then $(X\beta)'(X\beta_0)$$(X\beta)'(X\beta_0)=0$ for any $\beta_0 \in \text{nullspace}(C)$. This implies $X'X\beta \in \text{nullspace}(C)^\perp = \text{colspace}(C')$, so $X'X\beta=C'v$ for some $v$. Then, $X(X'X)^{-1} C' v = X\beta$, which shows $V_l \subseteq \text{colspace}(\tilde{X})$.
The more general case $\psi \ne 0$ can be obtained by slight modifications to the above proof. The fit of the restricted model would just be the projection $\tilde{P}_l$$\tilde{P}_r$ onto the affine space $\tilde{V}_r = \{X \beta : C \beta = \psi\}$, instead of the projection $P_r$ onto the subspace $V_r \{X\beta : C \beta = 0\}$$V_r =\{X\beta : C \beta = 0\}$. The two are quite related however, as one can write $\tilde{V}_l = V_r + \{X \beta_1\}$$\tilde{V}_r = V_r + \{X \beta_1\}$, where $\beta_1$ is an arbitrarily chosen vector satisfying $C \beta_1 = \psi$, and thus $$\tilde{P}_r y = P_r(y - X\beta_1) + X \beta_1.$$
Then, using the fact that $P_u X \beta_1 = X \beta_1$, we have $$(P_u - \tilde{P}_l) y = P_u y - P_r(y - X \beta_1) - X \beta_1 = (P_u - P_r)(y - X\beta_1) = P_l(y - X \beta_1).$$$$(P_u - \tilde{P}_r) y = P_u y - P_r(y - X \beta_1) - X \beta_1 = (P_u - P_r)(y - X\beta_1) = P_l(y - X \beta_1).$$ Recalling $P_l = \tilde{X} (\tilde{X}'\tilde{X})^{-1} \tilde{X}'$ with $\tilde{X} = X(X'X)^{-1} C'$, the difference in sum of squares of residuals can be shown to be $$(y - X \beta_1)' P_l (y - X \beta_1) = (C \hat{\beta} - \psi)'(C (X'X)^{-1} C') (C\hat{\beta} - \psi).$$