# Why do I find that OLS linear regression is robust against colinearity?

As per the textbook, OLS should fail when using colinear covariates. On their LinearRegression() documentation, sklearn states:

When features are correlated and the columns of the design matrix have an approximately linear dependence, the design matrix becomes close to singular and as a result, the least-squares estimate becomes highly sensitive to random errors in the observed target, producing a large variance. This situation of multicollinearity can arise, for example, when data are collected without an experimental design.

Why is it, then, that I have no problem doing a linear fit with perfectly colinear features? See below.

Note I am not talking about coefficients interpretations (into which colinearity obviously throws a monkey wrench), but simply the idea that colinearity + OLS leads to matrix singularity and error message.

import numpy as np
from scipy.stats import pearsonr
from sklearn.linear_model import LinearRegression
from matplotlib import pyplot as plt

x = (np.arange(0, 100, 1) + np.random.normal(1,1,100))
X = np.stack((x,x)).T

y = np.arange(0,200,2)

print(f"Example of X data: {X[:5]}")
print()
print(f"Example of y data: {y[:5]}")
print()

print(f"Correlation between the two features of X: {pearsonr(x, x)[0]}")
print()

model = LinearRegression()

model.fit(X, y)

fig, ax = plt.subplots(1, figsize=(10,5))

ax.scatter(x, y, label="truth")
ax.plot(x, model.predict(X), label="OLS")
ax.legend()
ax.set_xlabel("X")
ax.set_ylabel("y")

plt.show()

print(f"Regression intercept: {model.intercept_}")
print(f"Regression coefficients: {model.coef_}")


• There's no contradiction: multicollinearity translates to uncertainty in the coefficients, not in the fitted values. Your source unfortunately is ambiguous about this: by "the least-squares estimate" it means the coefficients only, not the fit.
– whuber
Commented Mar 23, 2022 at 13:39
• There IS a contradiction: by definition, OLS should not be able to produce a result with an input with perfectly colinear features, as OLS requires to compute the inverse of the Gram matrix of X, which is singular when X contains linearly dependent columns (and thus, non-invertible). Commented Mar 23, 2022 at 18:18
• You don't have perfectly correlated features.
– Dave
Commented Mar 23, 2022 at 18:28
• That's not a definition. I'm afraid it's an appeal to a couple of misconceptions. The definition is in the very name "ordinary least squares:" one minimizes the sum of squares of residuals. That problem always has a solution--but it cannot always be obtained by inverting a matrix. (The generalized inverse accomplishes this.) In the presence of collinearity, OLS has a space of solutions: but, by definition, they all minimize the sum of squared residuals and, because that objective is convex, the fit to the data is always the same for all solutions.
– whuber
Commented Mar 23, 2022 at 18:47

The issue with OLS is not that a solution is not achievable, but that the solution is not unique. For a model matrix, $$X$$, the quadratic form $$X^TX$$ needs to be invertible. The OLS coefficient is given by $$\hat{\beta} = \left(X^TX \right)^{-1}X^Ty$$. When $$X$$ is not of full rank, then a pseudo-inverse can be calculated instead. Whereas an inverse requires $$I = M^{-1}M$$, a pseudo-inverse is any matrix that has the property $$M = M^T M^{-1}M$$, but a pseudo-inverse is not unique. It looks like the particular expression in scipy is to "split" the beta between the two covariates. However, there are downstream issues of prediction, inference, and extrapolation that aren't as sensitive when the $$X$$ is not rank deficient.

• With the posted code X.shape returns (100,2) while y.shape returns (100,). This is a multivariate regression alright. Commented Mar 23, 2022 at 19:09
• Ah. I gave up Python for R before scipy came out, so I'm that rusty. I see the arange function for Y is a length 100 vector by 2s from 0 to 100. What scipy has done then is perform - and not tell you - a pseudoinverse. The solution is not unique and for some reason, it has split the coefficient value between the two X. Shame scipy doesn't warn you about this - R sets extra coefficients to NA - search for pseudoinverse in this SE to get more background. In other words, it's still a problem, but it doesn't mean a solution can't be found. Commented Mar 23, 2022 at 19:26
• Makes sense. I started looking into scipy's code and there seems to be a few addenda to basic OLS in case the input is not "proper." Thanks for pointing me in that direction AdamO! Commented Mar 23, 2022 at 19:44
• Subsequent question: non-colinearity amongst covariates is often stated to be one of OLS assumptions. Is it really? If so, why build an implementation of OLS that "works" (i.e. outputs a result) in that particular case? Commented Mar 24, 2022 at 14:35