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.