[proofplan]
We write the least squares problem through its normal equations, split into the $X$-block and the $Z$-block. The $Z$-block identifies the fitted contribution of $Z$ as the [orthogonal projection](/theorems/437) of $Y-X\hat{\beta}_X$ onto $\operatorname{Range}(Z)$, so applying $M_Z$ removes it. This reduces the $X$-block normal equation to $X^\top M_Z(Y-X\hat{\beta}_X)=0$. The full column rank of $[X\ Z]$ then gives invertibility of $X^\top M_ZX$, allowing us to solve for $\hat{\beta}_X$.
[/proofplan]
[step:Derive the block normal equations for the least squares minimizer]
Define the least squares objective function
\begin{align*}
Q: \mathbb{R}^k \times \mathbb{R}^m &\to \mathbb{R} \\
(\beta,\gamma) &\mapsto |Y - X\beta - Z\gamma|^2.
\end{align*}
Since $A=[X\ Z]$ has full column rank, the quadratic function $Q$ is strictly convex, and hence its minimizer $(\hat{\beta}_X,\hat{\gamma})$ is characterized by vanishing first derivatives in the $\beta$ and $\gamma$ directions.
For every $h \in \mathbb{R}^k$, differentiating $t \mapsto Q(\hat{\beta}_X+th,\hat{\gamma})$ at $t=0$ gives
\begin{align*}
0
&= -2h^\top X^\top(Y-X\hat{\beta}_X-Z\hat{\gamma}).
\end{align*}
Since this holds for all $h \in \mathbb{R}^k$,
\begin{align*}
X^\top(Y-X\hat{\beta}_X-Z\hat{\gamma})=0.
\end{align*}
Similarly, for every $\ell \in \mathbb{R}^m$, differentiating $t \mapsto Q(\hat{\beta}_X,\hat{\gamma}+t\ell)$ at $t=0$ gives
\begin{align*}
Z^\top(Y-X\hat{\beta}_X-Z\hat{\gamma})=0.
\end{align*}
Thus the block normal equations are
\begin{align*}
X^\top(Y-X\hat{\beta}_X-Z\hat{\gamma}) &= 0,\\
Z^\top(Y-X\hat{\beta}_X-Z\hat{\gamma}) &= 0.
\end{align*}
[/step]
[step:Use the $Z$-normal equation to eliminate $\hat{\gamma}$]
Because $[X\ Z]$ has full column rank, the columns of $Z$ are linearly independent. Hence $Z^\top Z$ is invertible. The second normal equation gives
\begin{align*}
Z^\top Y - Z^\top X\hat{\beta}_X - Z^\top Z\hat{\gamma}=0.
\end{align*}
Solving for $\hat{\gamma}$ yields
\begin{align*}
\hat{\gamma}
&= (Z^\top Z)^{-1}Z^\top(Y-X\hat{\beta}_X).
\end{align*}
Substituting this expression into the full residual gives
\begin{align*}
Y-X\hat{\beta}_X-Z\hat{\gamma}
&= Y-X\hat{\beta}_X
- Z(Z^\top Z)^{-1}Z^\top(Y-X\hat{\beta}_X)\\
&= M_Z(Y-X\hat{\beta}_X).
\end{align*}
[/step]
[step:Reduce the $X$-normal equation to the residualized regression equation]
Substitute the identity
\begin{align*}
Y-X\hat{\beta}_X-Z\hat{\gamma}=M_Z(Y-X\hat{\beta}_X)
\end{align*}
into the first normal equation. This gives
\begin{align*}
0
&= X^\top(Y-X\hat{\beta}_X-Z\hat{\gamma})\\
&= X^\top M_Z(Y-X\hat{\beta}_X)\\
&= X^\top M_ZY - X^\top M_ZX\hat{\beta}_X.
\end{align*}
Therefore
\begin{align*}
X^\top M_ZX\hat{\beta}_X = X^\top M_ZY.
\end{align*}
[/step]
[step:Prove that $X^\top M_ZX$ is invertible]
First note that $M_Z$ is symmetric and idempotent:
\begin{align*}
M_Z^\top &= M_Z,\\
M_Z^2 &= M_Z.
\end{align*}
Thus, for every $w \in \mathbb{R}^n$,
\begin{align*}
w^\top M_Zw
&= w^\top M_Z^\top M_Zw\\
&= |M_Zw|^2.
\end{align*}
We prove that $X^\top M_ZX$ has trivial kernel. Let $v \in \mathbb{R}^k$ satisfy
\begin{align*}
X^\top M_ZXv=0.
\end{align*}
Multiplying on the left by $v^\top$ gives
\begin{align*}
0
&= v^\top X^\top M_ZXv\\
&= |M_ZXv|^2.
\end{align*}
Hence $M_ZXv=0$. By the definition of $M_Z$,
\begin{align*}
Xv
&= Z(Z^\top Z)^{-1}Z^\top Xv.
\end{align*}
Define
\begin{align*}
a := (Z^\top Z)^{-1}Z^\top Xv \in \mathbb{R}^m.
\end{align*}
Then $Xv=Za$, so
\begin{align*}
[X\ Z]\begin{pmatrix}v\\-a\end{pmatrix}=0.
\end{align*}
Since $[X\ Z]$ has full column rank, its kernel is $\{0\}$. Therefore $v=0$ and $a=0$. Thus $X^\top M_ZX$ has trivial kernel. Since it is a $k \times k$ matrix, it is invertible.
[/step]
[step:Solve the residualized normal equation for the coefficient on $X$]
From the residualized normal equation,
\begin{align*}
X^\top M_ZX\hat{\beta}_X = X^\top M_ZY.
\end{align*}
Since $X^\top M_ZX$ is invertible, multiplying on the left by $(X^\top M_ZX)^{-1}$ gives
\begin{align*}
\hat{\beta}_X
&= (X^\top M_ZX)^{-1}X^\top M_ZY.
\end{align*}
This is exactly the coefficient vector obtained by regressing the residualized outcome $M_ZY$ on the residualized regressors $M_ZX$, completing the proof.
[/step]
