# scikit-learn linear regressor digests perfectly collinear features?

I am currently running this little piece of code:

import numpy as np
import sklearn.linear_model as skl_lm

X = np.random.normal(size = 1000)**2
y = X + np.random.normal(size = 1000)

Y = np.array([ [ X[i], X[i] ] for i in range(1000) ])

model = skl_lm.LinearRegression().fit(Y,y)
print(model.coef_)
print(np.dot(Y.transpose(),Y)/999)


It just creates 1000 squares of a normal random variable $$X$$, and uses it to create another linearly correlated variable $$y$$, accounting for some noise. Then creates a matrix of features $$Y$$, where the two features are actually both $$X$$. Then it tries to fit $$y$$ over $$Y$$. The model actually generates something, while the second print shows that the correlation matrix between columns of $$Y$$ is singular. How is sklearn.linear_model digesting this singular matrix? I believe it should not be able to determine the coefficients, as it should invert the covariance matrix.

Internally sklearn uses scipy.linalg.lstsq function for finding solutions to linear equation, which is the same as numpy.linalg.lstsq, and its innerworkings are described in this post: How does NumPy solve least squares for underdetermined systems?