Constrained Least Squares with Numpy and Python only

How to do constrained Least Squares only using numpy. Is there any way to incorporate constraints into numpy.linalg.lstsq() or is there any numpy+Python workaround to do constrained Least Squares? I know that I can do this easily with cvxpy and scipy.optimize, but this optimization will be a part of numba JITed function, and those two are not supported by numba.

EDIT:
Here is a dummy example of a problem I want to solve

arr_true = np.vstack([np.random.normal(0.3, 2, size=20),
                      np.random.normal(1.4, 2, size=20),
                      np.random.normal(4.2, 2, size=20)]).transpose()
x_true = np.array([0.3, 0.35, 0.35])
y = arr_true @ x_true

def problem():
    noisy_arr = np.vstack([np.random.normal(0.3, 2, size=20),
                        np.random.normal(1.4, 2, size=20),
                        np.random.normal(4.2, 2, size=20)]).transpose()

    x = cvxpy.Variable(3)
    objective = cvxpy.Minimize(cvxpy.sum_squares(noisy_arr @ x - y))
    constraints = [0 <= x, x <= 1, cvxpy.sum(x) == 1]
    prob = cvxpy.Problem(objective, constraints)
    result = prob.solve()
    output = prob.value
    args = x.value

    return output, args

Essentially, my problem is to minimize $Ax-y$ subject to $x_1+x_2+x_3 = 1$ and $\forall x_i, 0 \leq x_i \leq 1$.

Is it linear? Based on your numpy.linalg.lstsq suggestion I'll assume that for now. Also, I assume (to solve the simple case first) that you talk about equality constraints. Then you can easily solve it yourself by using Lagrange multipliers, and the solution can be found by solving a linear equation system.

If the problem is to minimize A * x - y over the constraints B * x = b, then solve

<img src="https://latex.codecogs.com/png.image?\large&space;\dpi{110}\bg{white}\begin{pmatrix}A^T\cdot&space;A&B^T\B&\end{pmatrix}\cdot\begin{pmatrix}x\\lambda\end{pmatrix}=\begin{pmatrix}A^T\cdot&space;y\b\end{pmatrix}" title="\begin{pmatrix}A^T\cdot A&B^T\B&\end{pmatrix}\cdot\begin{pmatrix}x\\lambda\end{pmatrix}=\begin{pmatrix}A^T\cdot y\b\end{pmatrix}" />

As an example:

from numpy import array, zeros
from numpy.linalg import solve

A = array([[1, 2], [3, 4], [5, 6], [7, 8]])
y = array([3.2, 6.8, 11.1, 15.2])

B = array([[1, -1]])
b = array([0.1])

S = zeros((3, 3))
r = zeros(3)

S[:2, :2] = A.T @ A
S[2:, :2] = B
S[:2, 2:] = B.T
r[:2] = A.T @ y
r[2] = b

x = solve(S, r)[:2]
print(x)

Disclamer: Bugs might be included, but brief testing gave me for x

[1.06262376 0.96262376]

This seems reasonable, and my constraint x1 - x2 = 0.1 is also fulfilled.