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?

  • $\begingroup$ Do you include a vector with 1,1...1 as first column in your X matrix? This is the intercept and is required to get a proper fit. $\endgroup$
    – Peter
    May 29 '19 at 15:52
  • $\begingroup$ Yes I do, in both cases. In fact in the 2D case the fit is pretty good $\endgroup$
    – Pigna
    May 29 '19 at 16:00
  • $\begingroup$ Your numbers don't look very big, but try to normalize your data (e.g. divide by the maximum) in any case. Linear regressions minimize the MSE ($(y-\hat y)^2$) therefore if you have a big y, say y=100, and a small y, say y=0, it will focus on the y=100 because the deviation has a higher impact on the residual. $\endgroup$ May 30 '19 at 17:28

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

from numpy.linalg import inv, det

print(det(X.T @ X))


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

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.