Old post but since I've been wondering about the same question here's my take on it. I came up a with the following numerical example to help me understand why the design matrix $\textbf{X}$ must be full rank.
First, with a linear model we want to solve (leaving out the error term $\boldsymbol{\epsilon}$):
$$\textbf{y} = \textbf{X} \boldsymbol{\beta}$$
where $\textbf{y}$ is vector of $n$ observation and $\textbf{X}$ is matrix of known covariates with $n$ rows and $p$ columns. We want solve to obtain the vector of $p$ coefficients $\boldsymbol{\beta}$.
Matrix $\textbf{X}$ is a transformation to pass from the vector space $\mathbb{R}^n$ of $\textbf{y}$ to the $\mathbb{R}^p$ space of $\boldsymbol{\beta}$. If $\textbf{X}$ is not full rank we cannot do this mapping uniquely, i.e. there is more than one solution for $\boldsymbol{\beta}$.
For example, with:
$$ \textbf{y} =
\begin{bmatrix}
8 \\
19 \\
27 \\
35 \\
\end{bmatrix}
\textbf{X} =
\begin{bmatrix}
1 & 2 \\
2 & 5 \\
3 & 7 \\
4 & 9 \\
\end{bmatrix}
$$
$\boldsymbol{X}$ is full rank and the model has a unique solution for $\boldsymbol{\beta} = \boldsymbol{(X^TX)^{-1}X^Ty} = [2 \ \ 3]$.
R code to reproduce this:
y <- c(8, 19, 27, 35)
X <- cbind(c(1, 2, 3, 4), c(2, 5, 7, 9))
qr(X)$rank == ncol(X) # Check X is full rank: TRUE
solve((t(X) %*% X)) %*% t(X) %*% y # Returns [2, 3]
# ---
# Alternatively, use the linear regression machinery:
lm(y ~ 0 + X[,1] + X[,2])
# Coefficients:
# X[, 1] X[, 2]
# 2 3
Now, consider:
$$ \textbf{y} =
\begin{bmatrix}
8 \\
16 \\
24 \\
32 \\
\end{bmatrix}
\textbf{X} =
\begin{bmatrix}
1 & 2 \\
2 & 4 \\
3 & 6 \\
4 & 8 \\
\end{bmatrix}
$$
$\boldsymbol{X}$ is not full rank (column two is twice column 1, i.e. it's linear combination of column 1) so there are multiple solutions for $\boldsymbol{\beta}$. For example with $\boldsymbol{\beta} = [2 \ \ 3]$, or $[4 \ \ 2]$, or $[6 \ \ 1]$. R code:
y <- c(8, 16, 24, 32)
X <- cbind(c(1, 2, 3, 4), c(2, 4, 6, 8))
qr(X)$rank == ncol(X) # FALSE: not full rank
X %*% c(2, 3) == y # All these solutions return TRUE
X %*% c(4, 2) == y
X %*% c(6, 1) == y
# Try to solve:
solve((t(X) %*% X)) %*% t(X) %*% y
Error in solve.default((t(X) %*% X)) :
Lapack routine dgesv: system is exactly singular: U[2,2] = 0