Use the standard results of partioned regression or equivalently FWL-theorem.
FWL-theorem
- The estimate of the coefficient $\beta_1$ in a regression model $$y = X_1 \beta_1 + X_2\beta_2 + e$$ will be the same as the estimate of $\beta_1$ in a regression model:
$$\bar y = \bar X_1 \beta_1 + e$$ where $\bar y$ are the residuals of regressing $y$ on $X_2$ and $\bar X_1$ are the residuals of regressing $X_1$ on $X_2$.
To do so start by considering the regression model
$$y = x\beta + e = X_1 \beta_1 + X_2\beta_2 + e$$
where define design matrix $X = [X_1 \ X_2]$ and $\beta = (\beta_1^\top \beta_2^\top)^\top$. The normal equations defining the OLS estimatior is given as
$$X^\top y = X^\top X \hat \beta $$ and written in partioned form this becomes
$$ \begin{bmatrix}X_1^\top y \\ X_2^\top y\end{bmatrix} = \begin{bmatrix}X_1^\top X_1 \hat \beta_1 + X_1^\top X_2\hat \beta_2 \\ X_2^\top X_1 \hat \beta_1 + X_2^\top X_2\hat \beta_2\end{bmatrix}$$
to find an expression for $\hat \beta_1$ use the lower part of this system of equations to get an expression for $\hat \beta_2$ which can be substituted into the upper part of the system. The lower part implies that
$$ (X_2^\top X_2)^{-1}X_2^\top y - (X_2^\top X_2)^{-1}X_2^\top X_1 \hat \beta_1 = \hat \beta_2$$ and inserting in the upper part results in the expression
$$ X_1^\top y = X_1^\top X_1 \hat \beta_1 + X_1^\top X_2[ (X_2^\top X_2)^{-1}X_2^\top y - (X_2^\top X_2)^{-1}X_2^\top X_1 \hat \beta_1 ]$$
defining the annihilator matrix (or residual maker matrix) $M_{X_2} :=I- X_2 (X_2^\top X_2)^{-1}X_2^\top$ this reduces to
$$ X_1^\top M_{X_2} y = X_1^\top M_{X_2}X_1 \hat \beta_1 $$
which can be solved to get an expression for $\hat \beta_1$
$$ ( X_1^\top M_{X_2}X_1 )^{-1}X_1^\top M_{X_2} y =\hat \beta_1 .$$
It is then noted that $M_{X_2}$ is idempotent and symmetric $M_{X_2} = M_{X_2}M_{X_2} = M_{X_2}^\top$ such that
$$ \hat \beta_1 = ( X_1^\top M_{X_2}X_1 )^{-1}X_1^\top M_{X_2} y = ( X_1^\top M_{X_2}^\top M_{X_2}X_1 )^{-1}X_1^\top M_{X_2}^\top M_{X_2} y$$
and then defining $\tilde X_1 := M_{X_2} X_1$ and $\tilde y := M_{X_2} y$ it follows that
$$ \hat \beta_1 = ( \tilde X_1^\top \tilde X_1 )^{-1}\tilde X_1^\top \tilde y$$
showing that $\hat\beta_1$ is the estimate you get if you regress $\tilde y$ on $\tilde X_1$. Furthermore the residual maker matrix $M_{X_2}$ is so called because it "computes" the residuals from a regression on $X_2$ so regressing $y$ on $X_2$ you get $\hat \lambda= (X_2^\top X_2)^{-1} X_2^{\top}y$ and predicted values $X_2\hat \lambda = X_2(X_2^\top X_2)^{-1} X_2^{\top}y$ and therefore resdiuals $y - X_2(X_2^\top X_2)^{-1} X_2^{\top}y = M_{X_2}y$. Similarly regressing the columns of $X_1$ on $X_2$ the same logic gives you residuals $M_{X_2}X_1$ which means that the estimate
$$ \hat \beta_1 = ( \tilde X_1^\top \tilde X_1 )^{-1}\tilde X_1^\top \tilde y$$
is the estimate of regressing residuals from a regression of $y$ on $X_2$ on the residuals of regression of the columns of $X_1$ on $X_2$.