Why does least squares linear regression perform so bad when switching from 2D to 3D line?

Question

The vector of coefficients that minimize least squares can be found like so:

beta = ((X'X)^-1)*X'y

Theoretically this result is true for any number of variables (columns of X).

I have made two python notebooks: one for univariate linear regression and another one for multivariate linear regression.

The code is analogous for both notebooks:

• create a random line

# multivariate snippet. A part from the plotting part,
# this is the only part of the code that differs between the two notebooks.
beta = np.random.randint(1,6,(3,))
X = np.linspace([1, 3, 7],[1, 14, 23],100)
y = X @ beta

• create noise

mu, sigma = 0, 1
noise = np.random.normal(mu, sigma, 100)

• add noise to y

y_noise = y + noise

• try to fit a line trough y_noise

beta_ = np.linalg.inv(X.T @ X) @ X.T @ y_noise
y_ = X @ beta_

• compare y_ with y_noise

Why do I get such a different result when fitting a 2D line and a 3D one?

Answer 1

The matrix $X^TX$ in your example is close to a singular:

from numpy.linalg import inv, det

print(det(X.T @ X))

1.479949989718308e-06

which makes its inverse very large

inv(X.T @ X)

[[ 4.83118746e+11  2.66548273e+11 -1.83251938e+11]
 [ 2.66548273e+11  1.47061116e+11 -1.01104517e+11]
 [-1.83251938e+11 -1.01104517e+11  6.95093558e+10]]

This is the reason of the error in your question. The following expression

inv(X.T @ X) @ X.T @ X

should produce an identity matrix. But in your example it is

[[ 1.2265625   1.10289171  4.8378709 ]
 [ 0.02282715  0.61713325  0.95782902]
 [-0.01519775  0.25675456  0.33339437]] 

So, when you multiply this matrix by beta the result is far from beta.

Another way to look at this singularity is that your 3D problem is not actually 3D. It is still 2D. That's why you have a line, not a 2D surface in your plot. It can be converted to a 2D representation by a transformation of coordinates. If you make a matrix for a proper 3D problem, say

np.random.seed(3)
X = np.hstack((np.ones((100,1)),np.random.rand(100,2)*100))

then inv(X.T @ X) @ X.T @ X gives the correct result:

[[ 1.00000000e+00 -1.24233956e-13 -6.21724894e-14]
 [-1.90819582e-17  1.00000000e+00  6.38378239e-16]
 [-8.67361738e-18  9.28077060e-16  1.00000000e+00]]

and the formula $\beta=(X^TX)^{-1}X^Ty$ works well.