No matrices, only scalars

Let’s, for the moment, forget about the matrices. Let’s say \(y = x * \beta\). We know \(x\) and \(y\) and we want to know \(\beta\). For this purpose, we multiply both sides of the equation by \(x^{-1}\). Put differently, then, \(x^{-1}*y = x^{-1}*x*\beta \implies x^{-1}*y=\beta \implies \beta = y/x\).

Special case: X as an invertible matrix

Before we move on to the more realistic, general, and common case where the feature matrix X can be non-invertible, let’s consider the case where X is invertible. The simplicity of the invertible matrices is that they, at least in the current case, can behave like scalars.

In this case, \(\mathbf{y}\) is a vector with, say, n elements, so y is a \(nx1\) matrix (a column vector with n elements) and \(\mathbf{X}\) is a \(nxn\) matrix (note that following the convention, we write matrices in bold). We still do not know \(\mathbf{\beta}\), and we know \(\mathbf{X}\) and \(\mathbf{y}\) (note that \(\mathbf{\beta}\) is also bold because for each feature, there is a different \(\beta\), so it has also become a vector). How do we find \(\mathbf{\beta}\) given \(\mathbf{y} = \mathbf{X} * \mathbf{\beta}\) ?

The first thing that would come to mind if we mistakenly thought that we are still in the real of scalars would be to divide \(\mathbf{y}\) by \(\mathbf{X}\). Unfortunately, that would not work, but we can use the matrix equivalent of the division operation; i.e. multiplying \(\mathbf{y}\) from the left by the inverse of \(\mathbf{X}\). So, we would get \(\mathbf{\beta} = \mathbf{X^{-1}} * \mathbf{y}\).

General case: X can be a non-invertible matrix

99.99% of the times, X would not be a square matrix. In fact, most of the times, X will have more observations (rows) than features/predictors (columns).¹ For a more intuitive grasp, let’s say we have collected sufficiently many observations, but we could only record a small number of features. Put concretely, let’s say that, \(\mathbf{y}\) has 100 determinants. If we knew all of these 100 determinants, we could always know the value of \(\mathbf{y}\) with perfect accuracy. But, in reality, we know only 5 of them. So, our feature matrix \(\mathbf{X}\) is \(100×5\) matrix instead of \(100×100\).

In this case, we can still try to guess what \(\mathbf{y}\) could be given the known 5 determinants. So, now, the formula became \(\mathbf{y} = \mathbf{X} * \mathbf{\beta} + \varepsilon\) . The last term is the error term. It is there because we do not know the other 95 determinants. But, for the moment, let’s ignore it and continue as if we still have \(\mathbf{y} = \mathbf{X} * \mathbf{\beta}\).

Because \(\mathbf{X}\) is not a square matrix, we cannot just take its inverse. However, what we can do is produce a square and invertible matrix using only \(\mathbf{X}\). The most straightforward way to do it is to multiply it from the left with its transpose. If \(\mathbf{X}\) is a \(100×5\) matrix, then \(\mathbf{X^{T}}*\mathbf{X}\) is a \(5×5\) matrix. So, we have \(\mathbf{X^{T}}*y = \mathbf{X^{T}}*\mathbf{X}*\beta\).

Still, \(\mathbf{\beta}\) is not alone, but now, we can take the inverse of what comes before \(\mathbf{\beta}\). Let’s do that:

\[(\mathbf{X^{T}}*\mathbf{X})^{-1}\mathbf{X^{T}}*y = (\mathbf{X^{T}}*\mathbf{X})^{-1}\mathbf{X^{T}}*\mathbf{X}*\beta\]

\[ \implies (\mathbf{X^{T}}*\mathbf{X})^{-1}\mathbf{X^{T}}*y = \beta\]

Voila!! We just arrived at the OLS estimator. Then, all this weird-looking formula is doing is trying to get a matrix that is invertible in front of \(\mathbf{\beta}\) and then leave \(\mathbf{\beta}\) alone.

The beauty of the formula is that it can be directly applied in the special cases as well:

If it is invertible, then this formula will produce \(\mathbf{X^{-1}}*y = \beta\), because individual elements in the parenthesis are also invertible. So, the formula, in the special case becomes \(\mathbf{X^{-1} * X^{T^{-1}} * X^{T} *y = \beta}\) and the transposed terms cancel each other out and yield \(\mathbf{X^{-1}*y = \beta}\).

If it is a scalar, because the transpose of a scalar is itself, the formula becomes \((1/x^{2}) * x * y = y/x = \beta\).

Let’s see these in action

Now, we will have two hats. When we wear the hat of the “Gods of Olympus”, we will know the underlying data generating process, read the truth. In fact, we will decide on the truth. When we wear the hat of the “researcher”, we will not know the truth, but we will see \(y\) or \(\mathbf{y}\), and \(x\) or \(\mathbf{x}\).

#Wearing the hat of the "Gods of Olympus"
set.seed(12345)
x <- rnorm(n=1)
true_beta <- 1
y <- x*true_beta

#Wearing the hat of the "researcher"

estimated_beta <- as.numeric(solve(t(x)%*%x) %*% t(x) %*% y)

all.equal(true_beta, estimated_beta)
## [1] TRUE

fitted_beta <- unname(lm(y~0+x)$coef)
all.equal(true_beta, fitted_beta)
## [1] TRUE

#Wearing the hat of the "Gods of Olympus"
X <- matrix(rnorm(100*100), nrow = 100)
true_beta <- matrix(seq(10,1, length.out = nrow(X)), ncol = 1)
y <- X %*% true_beta

#Wearing the hat of the "researcher"
estimated_beta <- solve(t(X)%*%X) %*% t(X) %*% y
all.equal(true_beta, estimated_beta)
## [1] TRUE

fitted_beta <- matrix(unname(lm(y~0+X)$coef), ncol = 1)
all.equal(true_beta, fitted_beta)
## [1] TRUE

observed_X <- X[,1:5]
estimated_noisy_beta <- solve(observed_X) %*% y
## Error in solve.default(observed_X): 'a' (100 x 5) must be square

(estimated_noisy_beta <- solve(t(observed_X) %*% observed_X) %*% t(observed_X) %*% y)
##           [,1]
## [1,] 13.373890
## [2,]  6.926159
## [3,]  3.317777
## [4,]  3.497302
## [5,]  7.713949

true_vs_est <- matrix(c(true_beta[1:5,], estimated_noisy_beta), byrow = FALSE, ncol = 2)
colnames(true_vs_est) <- c("truth", "estimated")
print(true_vs_est)
##          truth estimated
## [1,] 10.000000 13.373890
## [2,]  9.909091  6.926159
## [3,]  9.818182  3.317777
## [4,]  9.727273  3.497302
## [5,]  9.636364  7.713949