Since, there is a sense of ambiguity as to how the 'difference' of residuals of the restricted and unrestricted models leads to the $F$-statistic, it is deemed to be apt enough to briefly sketch the development in a general setting which would provide a better insight.
Theorem $[\rm I]:$ If $\mathbf x\sim \mathsf{MVN}(\boldsymbol\mu,\mathbf V),$
$$\mathbf x^\mathsf T\mathbf A\mathbf x\sim{\chi^2}^\prime\left[r(\mathbf A),\frac12\boldsymbol\mu^\mathsf T\mathbf A\boldsymbol \mu\right]\iff \mathbf{AV}~\text{idempotent}.^1$$
Consider the imposition of the constraint $$\mathbf K^\mathsf T\mathbf b = \mathbf m\tag 1$$ on the model $\mathbf y = \mathbf{Xb} +\mathbf e, $ where $\bf b$ is vector of parameters of order $k,~\mathbf K^\mathsf T $ is a matrix of order $s\times k$ along with the condition $r\left(\mathbf K^\mathsf T\right) = s, ~\mathbf m$ is a constant vector.
What's the effect of $(1) $ on the model? How should the estimator and the associated sum of squares look?
In oder to find the estimator of $\mathbf b,~\left(\tilde{\mathbf b}\right) $ subject to the constraint, one can employ Lagrange multiplier $(2\boldsymbol\lambda) $ in the minimisation of $$\left(\mathbf y -\mathbf X\tilde{\mathbf b}\right)^\mathsf T\left(\mathbf y -\mathbf X\tilde{\mathbf b}\right)+ 2\boldsymbol\lambda^\mathsf T\left(\mathbf K^\mathsf T\tilde{\mathbf b }-\mathbf m\right)\tag 2$$ w.r.t. $\tilde{\mathbf b}, ~\boldsymbol\lambda.$ Solving the equations
\begin{align}
\mathbf X^\mathsf T\mathbf X\tilde{\mathbf b} + \mathbf K\boldsymbol\lambda &= \mathbf X^\mathsf T\mathbf y\\ \mathbf K^\mathsf T\tilde{\mathbf b } &=\mathbf m,
\end{align}
yield
$$\tilde{\mathbf b }=\hat{\mathbf b }-(\mathbf X^\mathsf T\mathbf X)^{-1}\mathbf K\left[\mathbf K^\mathsf T(\mathbf X^\mathsf T\mathbf X)^{-1}\mathbf K\right]^{-1}\left(\mathbf K^\mathsf T\hat{\mathbf b }-\mathbf m\right).\tag 3$$
What would be the residual sum of squares? Computing $\left(\mathbf y -\mathbf X\tilde{\mathbf b}\right)^\mathsf T\left(\mathbf y -\mathbf X\tilde{\mathbf b}\right)$ using the fact that $\mathbf X^\mathsf T(\mathbf y -\mathbf X\hat{\mathbf b})= 0 $ would lead to
\begin{align}&= \left[\mathbf y -\mathbf X\hat{\mathbf b}+ \mathbf X\left(\hat{\mathbf b}-\tilde{\mathbf b}\right)\right]^\mathsf T\left[\mathbf y -\mathbf X\hat{\mathbf b}+ \mathbf X\left(\hat{\mathbf b}-\tilde{\mathbf b}\right)\right]\\&= \left( \mathbf y -\mathbf X\hat{\mathbf b}\right) ^\mathsf T \left( \mathbf y -\mathbf X\hat{\mathbf b}\right)+ \left(\hat{\mathbf b}-\tilde{\mathbf b}\right)^\mathsf T\mathbf X^\mathsf T\mathbf X\left(\hat{\mathbf b}-\tilde{\mathbf b}\right)\\ &\stackrel{(3)}{=} \textrm{SSE} +\underbrace{(\mathbf K^\mathsf T\hat{\mathbf b }-\mathbf m)^\mathsf T\left[\mathbf K^\mathsf T(\mathbf X^\mathsf T\mathbf X)^{-1}\mathbf K\right]^{-1}(\mathbf K^\mathsf T\hat{\mathbf b }-\mathbf m)}_{:= Q};\tag 4 \end{align}
so
$$\textrm{residual(reduced)} = \textrm{residual(full)} + Q. \tag 5^2$$
What is the distribution of $ Q? $
Notice that $\hat{\mathbf b}~\sim\mathsf{MVN}\left(\mathbf b, (\mathbf X^\prime\mathbf X) ^{-1}\sigma^2\right)$ as $\mathbf y \sim \mathsf{MVN}(\mathbf X\mathbf b,\sigma^2\mathbf I).$ Therefore $$\mathbf K^\mathsf T\hat{\mathbf b }-\mathbf m\sim\mathsf{MVN}\left(\mathbf K^\mathsf T\mathbf b-\mathbf m,\mathbf K^\mathsf T (\mathbf X^\prime\mathbf X) ^{-1} \mathbf K
\sigma^2\right).\tag 6$$ $Q$ can be seen to be a quadratic in $\mathbf K^\mathsf T\hat{\mathbf b}-\mathbf m, $ with $\left[\mathbf K^\mathsf T (\mathbf X^\prime\mathbf X) ^{-1} \mathbf K\right]^{-1}$ as the matrix of the quadratic. Now is the crucial part: apply Theorem $\rm[I] $ here (what are $\bf A, V$ here?) using $(6) $ to conclude that $$\frac Q{\sigma^2}\sim{\chi^2}^\prime\left[s,\frac{(\mathbf K^\mathsf T\hat{\mathbf b }-\mathbf m)^\mathsf T\left[\mathbf K^\mathsf T(\mathbf X^\mathsf T\mathbf X)^{-1}\mathbf K\right]^{-1}(\mathbf K^\mathsf T\hat{\mathbf b }-\mathbf m)}{2\sigma^2}\right].\tag 7$$ This is the desired result in a general framework. For constructing the required $F$-statistic, it has to be deduced that $Q$ and $\textrm{SSE}$ are independent.
Notes:
$[1] $ To prove this, compute the MGF of the quadratic form and use the property of the eigenvalues of an idempotent matrix.
$[2]$ It is tempting, perhaps, to interpret $Q,$ writing it in the form \begin{align}Q &= \mathbf y^\mathsf T\mathbf y -\text{SSE} -\left[ \mathbf y^\mathsf T\mathbf y - (\text{SSE}+ Q) \right]\\&=\text{SSR}- \left[ \mathbf y^\mathsf T\mathbf y - (\text{SSE}+ Q) \right]\\ &= \textrm{reduction(full)} - \left[ \mathbf y^\mathsf T\mathbf y - (\text{SSE}+ Q) \right],\tag{N1} \end{align} as the "reduction in sum of squares due to fitting of the reduced model" along the line of $(5);$ but this is not true, in general. In fact, the term in parenthesis in $\rm(N1) $ need not be even a sum of squares.
Reference:
Linear Models, S. R. Searle, John Wiley & Sons., $1971.$
