A = np.array([[1., 1.], [-1., -1.], [1., 0.]])

b = np.array([0., 0., 1.])

from scipy.sparse.linalg import lsqr
return_values = lsqr(A, b)
x = return_values[0]

A.T.dot(A).dot(x) == A.T.dot(b)

A [MxN] * x [Nx1] = b [Mx1]

U [MxM] * S[MxN] * V[N*N] * x[Nx1] = b[Mx1]

S[MxN] * V [N * N] * x[Nx1] = U.T [MxM] * b [Mx1]

>>> a = np.array([[1., 1.], [-1., -1.], [1., 0.]])
>>> b = np.array([0., 0., 1.])
>>> u, s, v = np.linalg.svd(a)
>>> np.allclose(u.T.dot(b)[-m+n:], 0) #check system is not overdetermined
True
>>> np.linalg.solve(s[:, None] * v, u.T.dot(b)[:n])
array([ 1., -1.])