Why is error of OLS not zero? Consider ordinary least squares (OLS): We have $n$ real datapoints $\mathbf{x}\in\mathbb{R}^d$ ($d$ features) organized in an ($n \times d$) matrix $X = (\mathbf{x}_1, \dots, \mathbf{x}_n)^T$ and we want to perform linear regression on $\mathbf{y} \in \mathbb{R}^n$ using a parameter vector $\boldsymbol{\theta} \in \mathbb{R}^d$. We write the loss function
\begin{equation}
L(\boldsymbol{\theta}) = \lVert \mathbf{y} - X\boldsymbol{\theta}\rVert^2 \tag{1}
\end{equation}
and set $\partial_\boldsymbol{\theta}L = 0$ to get the normal equations yielding the optimal parameter vector
\begin{equation}
\boldsymbol{\theta}^* = (X^T X)^{-1}X^T \mathbf{y} \tag{2}
\end{equation}
However, it seems to me that we can rearrange this using SVD: Write $X = U \Sigma V^T$ with orthogonal $n\times n$ matrix $U$, orthogonal $d\times d$ matrix $V$, and diagonal positive matrix $\Sigma$, then
\begin{align}
(X^T X)^{-1}X^T &= (V \Sigma^T U^T U \Sigma V^T)^{-1} V \Sigma^T U^T \tag{3a-e}
\\&= (V \Sigma^T \Sigma V^T)^{-1} V \Sigma^T U^T
\\&= V \Sigma^{-1} U^T
\\&= V \left(\Sigma^T U^T U (\Sigma^T)^{-1}\right)\Sigma^{-1}U^T
\\&= V \Sigma^T U^T (U \Sigma \Sigma^T U^T)^{-1}
\\&= X^T (X X^T)^{-1}
\end{align}
plugging (3e) into (2) gives $\boldsymbol{\theta}^* = X^T(X X^T)^{-1} \mathbf{y}$, and plugging this into (1) gives
\begin{align}
L(\boldsymbol{\theta}^*) = \lVert \mathbf{y} - X\left( X^T(X X^T)^{-1}\mathbf{y}\right)\rVert^2 = 0 \tag{4}
\end{align}
But this is nonsense, we can't possibly fit arbitrary $\mathbf{y}$ with a line. The only thing I can think of is that I'm somehow implicitly assuming $n < d$, but I can't see why this won't work for $n > d$.

My question is, why is this reasoning incorrect: what assumptions for OLS have I violated and what specific mathematical steps become invalid because of those errors?

How do you go from 3b to 3c?
You seem to be using the following $$(V\Sigma\Sigma V^T)^{-1}V\Sigma = V\Sigma^{-1}$$ But, that is only allowed when $\Sigma$ is square (ie when $n=p$). In that case $$(V\Sigma\Sigma V^T)^{-1} = (V^T)^{-1}\Sigma^{-1}\Sigma^{-1} V^{-1}$$
For $p=n$ it is correct that your cost function is zero.
