Roughly speaking: if in a linear regression the $p$ parameters are
subject to $m$ (independent) linear constraints, the model is simply a
reparameterisation of a linear model with $p- m$ unconstrained
parameters.
Consider the linear regression in standard matrix form
$$
\tag{1}
\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}.
$$
where $\mathbf{X}$ is $n \times p$. To stick to the OP notations, take $\beta_0 = \alpha$ and assume that the index
$i$ of $\beta_i$ runs from $0$ to $p-1$. The system of linear constraints
writes in matrix form as
$$
\tag{2} \mathbf{A}^\top \boldsymbol{\beta} = \mathbf{a}
$$
where $\mathbf{A}$ is a matrix with size $p \times m$ and $\mathbf{a}$ is a vector
of length $m$. We can assume
without loss of generality that $m < p$ and that $\mathbf{A}$ has rank
$m$. We can complete the $m$ columns of $\mathbf{A}$ with the $p-m$ columns
of a $p\times(p-m)$ matrix $\mathbf{B}$ so that the $p \times
p$ matrix
$$
\mathbf{C} = \begin{bmatrix}
\mathbf{A} &
\mathbf{B}
\end{bmatrix}
$$
is invertible. Then we can use $\boldsymbol{\gamma} := \mathbf{C}^\top
\boldsymbol{\beta}$ as a new parameter vector for the regression model
$\mathbf{Y} = \mathbf{Z}\boldsymbol{\gamma} +
\boldsymbol{\varepsilon}$ where $\mathbf{Z} :=
\mathbf{X}\mathbf{C}^{-\top}$, which is equivalent to (1). Using
blocks we can write $\mathbf{Z} =[\mathbf{Z}_1, \, \mathbf{Z}_2]$ where
the two blocks have $m$ and $p - m$ columns, and consistently
$\boldsymbol{\gamma}$ stacks $\boldsymbol{\gamma}_1$ and $\boldsymbol{\gamma}_2$. The
constraint (2) becomes $\boldsymbol{\gamma}_1 = \mathbf{a}$ in the new
parameterisation. So we just have to omit the $m$ first columns of
$\mathbf{Z}$ and fit the unconstrained linear regression model with parameter
$\boldsymbol{\gamma}_2$ and an offset
$$
\mathbf{Y} =
\underset{\text{offset}}{\underbrace{\mathbf{Z}_1 \mathbf{a}}} +
\mathbf{Z}_2\boldsymbol{\gamma}_2 +
\boldsymbol{\varepsilon},
$$
and then complete $\widehat{\boldsymbol{\gamma}}_2$ with
$\widehat{\boldsymbol{\gamma}}_1 = \mathbf{a}$ to get the full
$\widehat{\boldsymbol{\gamma}}$. The estimate
$\widehat{\boldsymbol{\beta}}$ for the original vector
$\boldsymbol{\beta}$ comes as $\widehat{\boldsymbol{\beta}} =
\mathbf{C}^{-\top} \widehat{\boldsymbol{\gamma}}$. We can infer on
$\boldsymbol{\gamma}_2$ as usual and derive predictions as well.
Now I shortly explain how to get $\mathbf{B}$ in practice. We can use
the QR decomposition of the $p \times m$ matrix $\mathbf{A}$
$$
\mathbf{A} = \underset{p \times p }{\underbrace{\mathbf{Q}}}
\overset{p \times m }{\overbrace{\mathbf{R}}}
=
\begin{bmatrix} \mathbf{Q}_1 \, \mathbf{Q}_2 \end{bmatrix}
\begin{bmatrix} \mathbf{R}_1 \\ \mathbf{0} \end{bmatrix} =
\mathbf{Q}_1 \mathbf{R}_1
$$
where $\mathbf{Q}$ is an orthogonal matrix and $\mathbf{R}$ is upper
triangular with non-zero diagonal elements. The blocks $\mathbf{Q}_1$
and $\mathbf{Q}_2$ have $m$ and $p-m$ columns, whith $\mathbf{R}_1$ is
$m \times m$ and invertible. We can take $\mathbf{B} :=
\mathbf{Q}_2$. Indeed, since $\mathbf{A} = \mathbf{Q}_1\mathbf{R}_1$,
the space $\text{cspan}(\mathbf{A})$ generated by the columns of
$\mathbf{A}$ in $\mathbb{R}^p$ is the same as
$\text{cspan}(\mathbf{Q}_1)$, with $\text{cspan}(\mathbf{Q}_2)$ being
a complement subspace.
For a R implementation, you need the design matrix
$\mathbf{X}$ and the matrix of restrictions $\mathbf{A}$. The qr.Q
function works on an object created by qr. It computes by default
only the thin QR decomposition $\mathbf{A} = \mathbf{Q}_1
\mathbf{R}_1$, and we have to use complete = TRUE to get the whole
$p \times p$ matrix $\mathbf{Q}$. Of course we could use the available
QR decomposition to avoid inverting $\mathbf{C}$ or its transpose.
## hypothetical design matrix
set.seed(135)
n <- 30; p <- 15; m <- 5
X <- cbind(1, matrix(runif(n * (p-1)), nrow = n))
## Constraints: mind the shift in the row indices due to the constant 'beta0'
A <- matrix(0.0, nrow = p, ncol = m)
a <- rep(0.0, m)
A[4, 1] <- A[5, 1] <- A[6, 1] <- 1.0
A[7, 2] <- A[8, 2] <- A[9, 2] <- 1.0
A[10, 3] <- A[13, 3] <- A[15, 3] <- 1.0
A[11, 4] <- A[13, 4] <- A[14, 4] <- 1.0
A[12, 5] <- A[14, 5] <- A[15, 5] <- 1.0
a[1] <- 1.0
## linear algebra
ind1 <- 1:m
ind2 <- (m + 1):p
QR <- qr(A)
Q <- qr.Q(QR, complete = TRUE)
B <- Q[ , ind2, drop = FALSE]
C <- cbind(A, B)
Cinv <- solve(C)
Z <- X %*% t(Cinv)
offset <- Z[ , ind1, drop = FALSE] %*% a
## Draw true coef and random observations
gamma2 <- rnorm(p - m)
Y <- Z %*% c(a, gamma2) + rnorm(n, sd = 0.1)
fit <- lm(Y ~ Z[ , ind2, drop = FALSE] - 1, offset = offset)
gammaHat <- c(a, coef(fit))
betaHat <- solve(t(C), gammaHat)
## check
t(A) %*% betaHat - a
```
