solving Ax =b for a non-square matrix A using python
Can use:

np.linalg.lstsq(x, y)
np.linalg.pinv(X) @ y
LinearRegression().fit(X, y)

Where 1 and 2 are from numpy and 3 is from sklearn.linear_model.
As a side note you will need to concatenate ones(use np.ones_like) in both 1 and 2 to represent the bias from the equation y = ax + b
I'm focusing on the special case where A is a n x d matrix (where k < d) representing an orthogonal basis for a subspace of R^d and b is known to be inside the subspace.
I thought of using the tools provided with numpy, however they only work with square matrices. I had the approach of filling in the matrix with some linearly independent vectors to "square" it and then solve, but I could not figure out how to choose those vectors so that they will be linearly independent of the basis vectors, plus I think it's not the only approach and I'm missing something that can make this easier.
is there indeed a simpler approach than the one I mentioned? if not, how indeed do I choose those vectors that would complete Ainto a square matrix?