The general form of a linear regression model is given by:
$$\hat{f}(x) = w_1 x_1 + \cdots + w_d x_d.$$Here $x_1, \ldots, x_d \in \mathrm{R}$ are input features, and $w_1, \ldots, w_d \in \mathrm{R}$ are model weights. In the real world, we'll be given a dataset of $N$ samples, and our job is to find the model weights that minimise the sum of squared errors. In particular, we want to minimise:
$$\sum_{i=1}^N \big( y_i - \hat{f}(x^{(i)}) \big)^2,$$where $y_1, \ldots, y_N \in \mathrm{R}$ are the targets, and $x^{(i)}$ represents the $i$-th sample. This can be expressed in matrix notation as:
$$\| y - Xw \|^2,$$where:
$$X = \begin{bmatrix} x^{(1)}_1 & x^{(1)}_2 & \cdots & x^{(1)}_d \\[3pt] x^{(2)}_1 & x^{(2)}_2 & \cdots & x^{(2)}_d \\[3pt] \vdots & \vdots & \ddots & \vdots \\[3pt] x^{(N)}_1 & x^{(N)}_2 & \cdots & x^{(N)}_d \\[3pt] \end{bmatrix} ,\quad w = \begin{bmatrix} w_1 \\[3pt] w_2 \\[3pt] \vdots \\[3pt] w_d \end{bmatrix} ,\quad y = \begin{bmatrix} y_1 \\[3pt] y_2 \\[3pt] \vdots \\[3pt] y_N \end{bmatrix}$$Our goal is to find a weight vector $w$ that minimises the least squares objective:
$$g(w) = \| y - Xw \|^2.$$We begin by noting the least squares objective function is convex. To see this, we find the Hessian of $g$ as:
$$\nabla^2 g(w) = X^TX,$$which is clearly positive semidefinite. Thus $g$ is a convex function. Convex differentiable functions have a special property - a point minimises the function if and only if the gradient at that point is zero. We find the gradient of $g$ as:
$$\nabla g(w) = 2X^TXw - 2X^Ty.$$Setting the gradient to zero, we obtain:
$$X^TXw = X^Ty.$$The above system of equations are called normal equations. Essentially, a point minimises $g$ if and only if it solves the normal equations
We’ve skipped most of the mathematical proofs and derivations in this section, as it’s meant to be brief. However, I’ll include some excellent references at the end of the article so you can dive deeper if you’re interested.
For illustrative purposes, let's assume the first two columns of $X$ are identical. So $X$ has linearly dependent columns. This implies $\operatorname{null} X \neq \{ 0 \}$. Intuitively, this means there exists a non-zero vector $v \in \mathrm{R}^d$ such that $Xv = 0$. In fact, In this case, we can find $v$. Define $v$ as follows:
$$v = \begin{bmatrix} -1 \\ 1 \\ 0 \\ \vdots \\ 0 \end{bmatrix},$$then we have:
$$\begin{align*} Xv &= -X_{:, 0} + X_{:, 1} + 0 \\[3pt] &= 0 \end{align*}$$Now, what does it imply when the null space of $X$ is non-empty? To explore this, let's assume $w$ satisfies the normal equations:
$$X^TXw = X^Ty.$$Then, for any $\alpha \in \mathrm{R}$, $w + \alpha v$ also satisfies the normal equations as illustrated below:
$$\begin{align*} X^TX(w + \alpha v) &= X^TXw + \alpha X^TXv \\[4pt] &= X^TXw + \alpha X^T(Xv) \\[4pt] &= X^TXw + \alpha X^T(0) \\[4pt] &= X^TXw + 0 \\[4pt] &= X^TXw \end{align*}$$On a side note, we need not assume that a solution to the normal equations exists. It can be mathematically proven that a solution $w$ is always guaranteed, and a formula for it can be derived using the Singular Value Decomposition (SVD) of $X$. While we will not delve into the details here, I'll provide a resource for the interested reader at the end of this article.
We've just seen that the normal equations can be satisfied by infinitely many vectors, which in turn implies that the objective function
$$g(w) = \| y - Xw \|^2$$can be minimised by infinitely many vectors. In essence, the solution isn't unique.
Suppose you've found two weights $w_1$ and $w_1 + \alpha v$ that satisfies the normal equations and want to deploy one of them to production. Does it matter which one you deploy? Before you deploy, you want to measure their individual performance on a test set $X^\prime$. Here's a crucial point: If the first two columns of the test set are identical as well, then both $w_1$ and $w_1 + \alpha v$ will give the same results. However, if the columns are not identical, as it often happens in the real world, then the results will be way off.
So, in essence, if both the train set features and test set features are correlated in the same way, then multicollinearity wouldn't be a problem. If they're not, then we'd have to deal with an unstable and uninterpretable model.
The model is uninterpretable because there are infinitely many weights to choose from. Each of these weights tells a different story about feature contribution, yet they all minimize the objective function.
Since infinitely many solutions exist, it's possible that we have a picked a model with huge weights. This happens if $\alpha$ is large. In general, a model with huge weights is unstable because a small change in the values of the features will lead to wildly different outcomes. So it's always a good idea to pick a model with small weights. This can be done through regularisation