<s>[INST] Why is $n < p$ a problem for OLS regression? I realize I can't invert the $X'X$ matrix but I can use gradient descent on the quadratic loss function and get a solution. I can then use those estimates to calculate standard errors and residuals. Am I going to encounter any problems doing this?
 [/INST] Here is a little specific example to illustrate the issue:
Suppose you want to fit a regression of $y_i$ on $x_i$, $x_i^2$ and a constant, i.e. 
$$ y_i = a x_i + b x_i^2 + c + u_i $$
or, in matrix notation,
\begin{align*}
 \mathbf{y} =
 \begin{pmatrix}
  y_1 \\
  \vdots \\
  y_n 
 \end{pmatrix}, \quad
 \mathbf{X} =
 \begin{pmatrix}
  1 & x_1 & x_1^2 \\
  \vdots & \vdots & \vdots \\
  1 & x_n & x_n^2
 \end{pmatrix}, \quad
 \boldsymbol{\beta} =
 \begin{pmatrix}
  c \\ a \\ b
 \end{pmatrix}, \quad
 \mathbf{u} =
 \begin{pmatrix}
  u_1 \\ \vdots \\ u_n
 \end{pmatrix}
\end{align*}
Suppose you observe $\mathbf{y}^T=(0,1)$ and $\mathbf{x}^T=(0,1)$, i.e., $n=2<p=3$. 
Then, the OLS estimator is
\begin{align*}
 \widehat{\boldsymbol{\beta}} =& \, (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y} \\
 =& \, \left[
    \begin{pmatrix}
   1 & 0 & 0 \\
   1 & 1 & 1  
    \end{pmatrix}^T
    \begin{pmatrix}
     1 & 0 & 0 \\
   1 & 1 & 1
    \end{pmatrix}
    \right]^{-1}
    \begin{pmatrix}
   1 & 0 & 0 \\
   1 & 1 & 1  
    \end{pmatrix}^T
    \begin{pmatrix}
     0 \\ 1
    \end{pmatrix} \\
 =& \, \left[\underbrace{\begin{pmatrix}
   2 & 1 & 1 \\
   1 & 1 & 1 \\
   1 & 1 & 1
    \end{pmatrix}}_{\text{not invertible, as $\mathrm{rk}()=2\neq 3$}}
\right]^{-1}  
    \begin{pmatrix}
     1 \\ 1 \\ 1
    \end{pmatrix}
\end{align*}
There are infinitely many solutions to the problem: setting up a system of equations and inserting the observations gives
\begin{align*}
 0 =& \, a \cdot 0 + b \cdot 0^2 + c \ \ \Rightarrow c=0 \\
 1 =& \, a \cdot 1 + b \cdot 1^2 + c \ \ \Rightarrow 1 = a + b
\end{align*}
Hence, all $a=1-b$, $b \in \mathbb{R}$ satisfy the regression equation.
</s>