Why does this multiple linear regression fail to recover the true coefficients?

Question

I am trying to use linear least squares regression to extract the coefficients of a model. Specifically, I am looking at a model with two independent predictor variables $x_1$ and $x_2$, and an output response variable $y$, with coefficients $\beta_0,\beta_1$ and $\beta_2$. (I believe this is a case of multiple linear regression.) The model is of the form $$ y_i = \beta_0 + \beta_1x_{i1} + \beta_2 x_{i2} + \varepsilon_i \tag{1} $$ where $i$ denotes the observation number. Or in matrix form with $N$ observations $$ \begin{align} \mathbf{Y} &= \mathbf{A}\boldsymbol{\beta} + \boldsymbol{\varepsilon} \\ \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_N \end{pmatrix} &= \begin{pmatrix} 1 & x_{11} & x_{21} \\ 1 & x_{12} & x_{22}\\ \vdots & \vdots & \vdots \\ 1 & x_{1N} & x_{2N} \end{pmatrix} \begin{pmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \end{pmatrix} + \begin{pmatrix} \varepsilon_1 \\ \varepsilon_2\\ \vdots \\ \varepsilon_N \end{pmatrix} \end{align} \tag{2} $$ I should be able to solve this using the least squares condition $$\boldsymbol{\beta} = (\mathbf{A}^\textrm{T}\mathbf{A})^{-1}\mathbf{A}^\textrm{T}\mathbf{Y} \tag{3}$$.

I made a simple (Matlab) script to see if I could recover the correct beta coefficients on a test case, given below. I used three methods: 1) directly solving Eq.(3) using Matlab's backslash operator 2) the same thing, except neglecting the transpose in Eq.(3), since it coincides with the least squares result. 3) solving Eq.(3) using a QR decomposition insead.

x1_data = linspace(-1,1,20).';      % First independent variable
x2_data = linspace(-0.3,0.6,20).';  % Second independent variable

beta0 = -0.2; % Intercept
beta1 = 0.4;  % coefficient of x1 data
beta2 = 1.2;  % coefficient of x2 data

y_data =  beta0 + beta1*x1_data + beta2*x2_data; % Create test response data

A = [ ones(length(x1_data),1)   x1_data   x2_data ]; % Design matrix
Y = y_data;

beta_fit1 = (A.'*A) \ (A.'*Y);  % Method #1
beta_fit2 = A \ Y;              % Method #2

[Q,R] = qr(A);
beta_fit3 = R \ (Q'*Y);         % Method #3

which outputs:

beta_fit1 = [-254.4 -762.3 1696.0]
beta_fit2 = [-0.02 0.94 0]
beta_fit= [-0.02 0.94 0]

None of which match the original input beta coefficients. Can anyone tell me what might be going wrong here?

Thanks

---

Edit

An example to show that it works when fitting a 2D parabolic function with no noise (i.e. with deterministic data, as mentioned by @Gregg H):

x1_data = linspace(-1,1,20).';      % First predictor variable
x2_data = linspace(-0.3,0.6,20).';  % Second predictor variable

beta0_true = -0.2; % Intercept
beta1_true = 0.4;  % coefficient of x1 data
beta2_true = 1.2;  % coefficient of x2 data

y_data =  beta0_true + beta1_true*x1_data.^2 + beta2_true*x2_data.^2 ; % Create test response data

A = [ ones(length(x1_data),1)   x1_data.^2   x2_data.^2 ]; % Design matrix
Y = y_data;

beta_fit2 = A \ Y;

which outputs as expected:

beta_fit = [-0.2   0.4   1.2]

The issue is when trying to add a linear term (such as beta4_true * x_data1 ), i.e. a tilt to the parabola.

Answer 1

At least the 2nd and 3rd solutions are correct.

Your design matrix has dependent variables. For example the third column can be expressed in terms of the first two columns $x_3 = 0.15 + 0.45 x_2$ and the equation can also be expressed as

$$\begin{array}{} -0.2 + 0.4 x_ 1 + 1.2 x_3 &=& -0.2 + 0.4 x_ 1 + 1.2 (0.15 + 0.45 x_2) \\ &=& -0.2 + 0.4 x_ 1 + 1.2 (0.15 + 0.45 x_2)\\ &=& -0.02 + 0.94 x_1 \end{array}$$

Methods 2 and 3

This last equation on the right hand side is the solution given by the 2nd and 3rd methods which probably drop one of the columns. In R you get the same behavior when we use the function lm which gives as output

> lm(y~X+0)

Call:
lm(formula = y ~ X + 0)

Coefficients:
   X1     X2     X3  
-0.02   0.94     NA  

The last column is ignored when you give a computer the task to solve the equation.

Method 1

Your 1st method probably attempts to inverse the (non-invertible) matrix anyway. For example, the inverse command does give some output. In my computer (an online https://www.tutorialspoint.com/execute_matlab_online.php) I get:

disp(inv(A.'*A));
-3.3777e+13  -1.0133e+14   2.2518e+14
-1.0133e+14  -3.0399e+14   6.7554e+14
 2.2518e+14   6.7554e+14  -1.5012e+15

and gives some output that is close but not exact (possibly due to round of errors).

In my case I got -0.071429 0.785714 0.342857, which is close to a correct solution $-0.071429+0.15 \cdot 0.342857 \approx -0.02$ and $0.785714+0.45 \cdot 0.342857 \approx 0.94$

In your case the difference is larger $-254.4 + 1696\cdot 0.15 \approx 0$ and $-762.3 + 1696*0.45 = 0.9$ (but this might be due to the output being given with less precision)

In R I can get the same result when I use the solve command while setting the tolerance parameter extremely low. In that case the inverse matrix is still computed; and it can be computed because the columns in the matrix X are not entirely dependent due to round off errors.

X = cbind(rep(1,20), seq(-1,1,length.out = 20), seq(-0.3,0.6,length.out = 20))
beta = c(-0.2,0.4,1.2)
y = X %*% beta
X = round(X,5)
solve(t(X) %*% X, tol = 10^-50) %*% t(X) %*% y
#           [,1]
#[1,] -0.0430674
#[2,]  0.8707996
#[3,]  0.1537806